#     This file is part of the Test Set for BVP solvers
#
#        http://www.dm.uniba.it/~bvpsolvers/
#
#        Problem bvpT6
#        ODE of dimension 2
#
#     DISCLAIMER: see
#     http://www.dm.uniba.it/~bvpsolvers/testsetbvpsolvers
#
#     
# this example require the library pracma  for the function erf()


bvpT6<- function(){

  
  
prob<- function(){ 
       fullnm <- 'Problem bvpT6'
       problm = 'bvpT6'
       typebvp   = 'SPBVP'
    
       neqn = 2
       nlbc = 1
       aleft    = -1.0e0
       aright   = 1.0e0
       ms <- rep(1,neqn)
       numjac   = FALSE
       numbcjac = FALSE
       linear   = TRUE
       Rpar<-c(0.5)
       Ipar<-0
       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))
}

# pracma is a library for use the function erf()

library(pracma)
  
  

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))
}
  
  

feval = function(x,y,eps,Rpar,Ipar){ 
        f <- c(y[2],
               (-x*y[2]-eps*pi*pi*cos(pi*x)-pi*x*sin(pi*x))/eps )
    
        return(list(f))
}
  


jeval = function(x,y,eps,Rpar,Ipar){
        dfy <- matrix( nrow=2,ncol=2,byrow = TRUE, data = c(
                       0, 1,
                       0, -x/eps))
    
        return((dfy))             
    
}
  
 
bceval <- function(i, y, eps,Rpar,Ipar) {
    
         if (i == 1) return(y[1]+2)
         if (i == 2) return(y[1])
    
    
}
  
  
dbceval <- function(i, y, eps,Rpar,Ipar) {
    
          if (i == 1) return(c(1, 0))
          if (i == 2) return(c(1, 0))
    
}
  
  
  
setoutput<- function(neqn,plotsol=NULL){ 
            solref = TRUE
    
            if (is.null(plotsol)){ 
                 nindsol = neqn
                 indsol = 1;neqn}
            else{
                 nindsol = length(plotsol)
                 indsol = plotsol}         
            return(list(solref=solref,nindsol=nindsol,indsol=indsol))
} 
  


esolu <- function(X,parms,Rpar,Ipar){
         lambda=parms
         Sqep<-sqrt(2.e0*lambda);
    
         Exact =matrix( c( cos(pi*X)+erf(X/Sqep)/erf(1.e0/Sqep), 
                           2^(1/2)/(pi^(1/2)*lambda^(1/2)*exp(X^2/(2*lambda))*erf(2^(1/2)/(2*lambda^(1/2)))) - pi*sin(pi*X)),
                           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))   
  
}