#     This file is part of the Test Set for BVP solvers
#
#        http://www.dm.uniba.it/~bvpsolvers/
#
#        Problem bvpT7
#        BVP of dimension 2
#
#     DISCLAIMER: see
#     http://www.dm.uniba.it/~bvpsolvers/testsetbvpsolvers
#
#
bvpT7<- function(){

  
prob<- function(){ 
       library(pracma)
       fullnm <- 'Problem bvpT7'
       problm = 'bvpT7'
       typebvp   = 'SPBVP'
    
       neqn = 2
       nlbc = 1
       aleft    = -1.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))
}
  
  
  
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]+y[1]-(1.0e0+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,
                       1/eps, -x/eps))
    
        return((dfy))             
    
}
  


bceval <- function(i, y, eps,Rpar,Ipar) {
      
         if (i == 1) return(y[1]+1)
         if (i == 2) return(y[1]-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
   
              Aa <-X*X/(2.0e0*lambda)
              Bb <- 1.0e0/(2.0e0*lambda)
              Cc <- exp(-Aa)
              Dd <- exp(-Bb)

    
              Exact =matrix( c( cos(pi*X)+X+(X*erf(X/sqrt(2.0e0*lambda))+sqrt(2.0e0*lambda/pi)*Cc)/
                            (erf(1.0e0/sqrt(2.0e0*lambda))+sqrt(2.0e0*lambda/pi)*Dd),
                             erf((2^(1/2)*X)/(2*lambda^(1/2)))/(erf(2^(1/2)/(2*lambda^(1/2))) + (2^(1/2)*lambda^(1/2))/(pi^(1/2)*exp(1/(2*lambda)))) - pi*sin(pi*X) + 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))   
  
}