#     This file is part of the Test Set for BVP solvers
#
#        http://www.dm.uniba.it/~bvpsolvers/
#
#        Problem bvpT1
#        SPBVP of dimension 2
#
#     DISCLAIMER: see
#     http://www.dm.uniba.it/~bvpsolvers/testsetbvpsolvers/?page_id=34
#
#     The most recent version of this source file can be found at
#     http://www.dm.uniba.it/~bvpsolvers/bvpTestSet/Rsrc/problems/bvpT1.R
#
#
#


bvpT1 <- function(){
  
#------ function prob -------#  
prob<- function(){ 
    fullnm <- 'Problem bvpT1'
    problm  = 'bvpT1'
    typebvp = 'SPBVP'
    neqn = 2
    nlbc = 1
    aleft    = 0.0e0
    aright   = 1.0e0
    numjac   = FALSE
    numbcjac = FALSE
    linear   = TRUE
    Rpar<-c(0.5)
    Ipar<-0
    ms <- rep(1,neqn)
   
      
    return(list(fullnm=fullnm,problm=problm,typebvp=typebvp,neqn=neqn,
                 nlbc=nlbc,ms=ms,aleft=aleft,aright=aright,Rpar=Rpar,Ipar=Ipar,
                 numjac=numjac,numbcjac=numbcjac,linear=linear))
}
  
#----- function init -----#  
  
init <- function(neqn,ms,aleft,aright) { 
    givmsh = FALSE
    givey = FALSE
    nmsh = 11
    xguess = seq(aleft,aright,by=(aright-aleft)/(nmsh-1))
    
    yguess <- NULL
    for (i in 1:neqn){ 
       for (j in 1:ms[i])
          yguess <- rbind(yguess,rep(0,nmsh))
    }
    
  
    return(list(givmsh=givmsh,givey=givey,
                nmsh=nmsh,xguess=xguess,yguess=yguess))
}
  

#----- function feval -----#  

feval = function(x,y,eps,Rpar,Ipar){ 
           f <- c(y[2],
                  y[1]/eps)
        return(list(f))
}
  
  
#----- function jeval -----#  

jeval = function(x,y,eps,Rpar,Ipar){
                dfy <- matrix( nrow=2,ncol=2,byrow = TRUE, data = c(
                0, 1,
                1/eps, 0))
        return((dfy))             
}
  
  
#----- function bceval -----#  

bceval <- function(i, y, eps,Rpar,Ipar) {
  
                  if (i == 1) return(y[1]-1)
                  if (i == 2) return(y[1])
   
}
  

#----- function dbceval -----#  

dbceval <- function(i, y, eps,Rpar,Ipar) {
    
                   if (i == 1) return(c(1, 0))
                   if (i == 2) return(c(1, 0))
    
}
  
  
#----- function  setoutput -----#  
  
setoutput<- function(neqn,plotsol=NULL){ 
                     solref = TRUE
                     if (is.null(plotsol)){ 
                         nindsol = 1
                         indsol = 1}
                    else{
                         nindsol = length(plotsol)
                         indsol = plotsol}         
                    return(list(solref=solref,nindsol=nindsol,indsol=indsol))
}
  
  
#----- function esolu -----#  

esolu <- function(X,parms,Rpar,Ipar){
                  lambda=parms
                  print(lambda)
                  Exact =matrix( c( (exp(-X/sqrt(lambda))-exp(-(2.e0-X)/sqrt(lambda)))/(1.e0-exp(-2.e0/sqrt(lambda))),
                  (1/(lambda^(1/2)*exp(X/lambda^(1/2))) + exp((X - 2)/lambda^(1/2))/lambda^(1/2))/(1/exp(2/lambda^(1/2)) - 1)), 
                  ncol=2,nrow=length(X),byrow=FALSE)
  
                  return(Exact=Exact)}
  
return(list(prob=prob,init=init,feval=feval,jeval=jeval,bceval=bceval,dbceval=dbceval,setoutput=setoutput,esolu=esolu))   
  
}