#     This file is part of the Test Set for BVP solvers
#
#        http://www.dm.uniba.it/~bvpsolvers/
#
#        Problem bvpT20
#        ODE of dimension 2
#
#     DISCLAIMER: see
#     http://www.dm.uniba.it/~bvpsolvers/testsetbvpsolvers
#


bvpT20<- function(){
  
  
prob<- function(){ 
        fullnm <- 'Problem bvpT20'
        problm = 'bvpT20'
        typebvp   = 'SPBVP'
        
        neqn = 2
        nlbc = 1
        aleft    = 0
        aright   = 1
        numjac   = FALSE
        numbcjac = FALSE
        linear   = FALSE
        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],
              (1.e0-y[2]*y[2])/eps)
    
        return(list(f))
}
  
  
jeval = function(x,y,eps,Rpar,Ipar){
        dfy <- matrix( nrow=2,ncol=2,byrow = TRUE, data = c(
                        0, 1.0e0,
                        0,-2.e0*y[2]/eps))
    
        return((dfy))             
}
  
  

bceval <- function(i, y, eps,Rpar,Ipar) {

          if (i == 1) {x =-0.745e0/eps
                          return (y[1]-1.e0-eps*(-x+log((exp(2.e0*x)+1.e0)/2.e0)))
                        }
          if (i == 2) {X = 0.255e0/eps
                          return( y[1]-1.e0-eps*(X+log((exp(-2.e0*X)+1.e0)/2.e0)))
                        }
}
  
  
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))
}
  

# todo the exact solution should be corrected
 
esolu <- function(X,parms,Rpar,Ipar){
          lambda=parms
          Xx<-(X-0.745e0)/lambda 
          Exact =matrix( c(0,0),
                         ncol=2,nrow=length(X),byrow=FALSE)
          for (j in 1:length(X))
            if(Xx[j]>0){
              Exact[j,] =matrix( c(1.e0+lambda*(Xx[j]+log((1.e0+exp(-2.e0*Xx[j]))/2.e0)),
                           lambda*(1/lambda - 1/(lambda*exp((2*(X[j] - 149/200))/lambda)*(1/(2*exp((2*(X[j] - 149/200))/lambda)) + 1/2)))),
                           ncol=2,nrow=1,byrow=FALSE)
                    }
            else{
              Exact[j,] =matrix( c( 1.e0+lambda*(-Xx[j]+log((1.e0+exp(2.e0*Xx[j]))/2.e0)),
                           -lambda*(1/lambda - exp((2*(X[j] - 149/200))/lambda)/(lambda*(exp((2*(X[j] - 149/200))/lambda)/2 + 1/2)))),
                            ncol=2,nrow=1,byrow=FALSE)
                    }
      
    
          return(Exact=Exact)
      }
  
  return(list(prob=prob,init=init,feval=feval,jeval=jeval,bceval=bceval,dbceval=dbceval,setoutput=setoutput,esolu=esolu))   
  
}