#     This file is part of the Test Set for BVP solvers
#
#        http://www.dm.uniba.it/~bvpsolvers/
#
#        Problem bvpT31
#        ODE of dimension 4
#
#     DISCLAIMER: see
#     http://www.dm.uniba.it/~bvpsolvers/testsetbvpsolvers
#
#       Authors:
#
#       Jeff R. Cash 
#            (Department of Mathematics, Imperial College,  London, England.)
#       Davy  Hollevoet 
#            (Vakgroep Toegepaste Wiskunde en Informatica, Universiteit Gent, Belgium.)
#       Francesca Mazzia  
#            (Dipartimento di Matematica, Universita' di Bari, Italy)
#       Abdelhameed Nagy Abdo
#            (Dipartimento di Matematica, Universita' di Bari, Italy)
#            (Dept. of Mathematics, Faculty of Sciences, Benha  University,Egypt)


bvpT31<- function(){
 
  
prob<- function(){ 
        fullnm <- 'Problem bvpT31'
        problm = 'bvpT31'
        typebvp   = 'SPBVP'
        
        neqn = 4
        nlbc = 2
        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(sin(y[2]),
               y[3],
               -y[4]/eps,
               ((y[1]-1.e0)*cos(y[2])-y[3]*(1.e0/cos(y[2])+eps*y[4]*tan(y[2])))/eps)
    
        return(list(f))
}
  

jeval = function(x,y,eps,Rpar,Ipar){
    
        Z2<-cos(y[2])
        Z3sc<- 1.e0/Z2
        dfy <- matrix( nrow=4,ncol=4,byrow = TRUE, data = c(
                       0, cos(y[2]), 0, 0,
                       0,0,1.0e0 , 0,
                       0, 0, 0,  -1.0e0/eps,
                       Z2/eps,(-(y[1]-1.e0)*sin(y[2])-y[3]*Z3sc*(tan(y[2])+eps*y[4]*Z3sc))/eps, -(Z3sc+eps*y[4]*tan(y[2]))/eps, (-y[3]*eps*tan(y[2]))/eps ))
          
        return((dfy))             
    
}
  

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

          if (i == 1) return(y[1])
          if (i == 2) return(y[3])
          if (i == 3) return(y[1])
          if (i == 4) return(y[3])  
}
  
  

dbceval <- function(i, y, eps,Rpar,Ipar) {
    
          if (i == 1) return(c(1, 0, 0, 0))
          if (i == 2) return(c(0, 0, 1, 0))
          if (i == 3) return(c(1, 0, 0, 0))
          if (i == 4) return(c(0, 0, 1, 0))
    
}
  
  
setoutput<- function(neqn,plotsol=NULL){ 
              solref = FALSE
              
              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
          Exact =matrix( c(),
                         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))   
  
}