#     This file is part of the Test Set for BVP solvers
#
#        http://www.dm.uniba.it/~bvpsolvers/
#
#        Problem bvpT33
#        ODE of dimension 6
#
#     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)


bvpT33<- function(){

  
  
prob<- function(){ 
        fullnm <- 'Problem bvpT33'
        problm = 'bvpT33'
        typebvp   = 'SPBVP'
        
        neqn = 6
        nlbc = 3
        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 = TRUE
        givey = TRUE
        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.5,nmsh))
        return(list(givmsh=givmsh,givey=givey,
                        nmsh=nmsh,xguess=xguess,yguess=yguess))
}
  
feval = function(x,y,eps,Rpar,Ipar){ 
        f <- c(y[2],
               (y[1]*y[4]-y[3]*y[2])/eps,
               y[4],
               y[5],
               y[6],
               (-y[3]*y[6]-y[1]*y[2])/eps)
           
    
        return(list(f))
}
  
jeval = function(x,y,eps,Rpar,Ipar){
        dfy <- matrix( nrow=6,ncol=6,byrow = TRUE, data = c(
                        0, 1.0e0, 0, 0, 0, 0,
                        y[4]/eps, -y[3]/eps , -y[2]/eps, y[1]/eps, 0, 0,
                        0, 0, 0, 1.0e0, 0, 0,
                        0, 0, 0, 0, 1.0e0, 0,
                        0, 0, 0, 0, 0, 1.0e0,
                        -y[2]/eps, -y[1]/eps, -y[6]/eps, 0, 0, -y[3]/eps ))
    
        return((dfy))             
}
  
  bceval <- function(i, y, eps,Rpar,Ipar) {
    
            if (i == 1) return(y[1]+1.0e0)
            if (i == 2) return(y[3])
            if (i == 3) return(y[4])
            if (i == 4) return(y[1]-1-0e0)
            if (i == 5) return(y[3])
            if (i == 6) return(y[4])
}
  
  

dbceval <- function(i, y, eps,Rpar,Ipar) {
    
            if (i == 1) return(c(1, 0,0,0,0,0))
            if (i == 2) return(c(0, 0,1,0,0,0))
            if (i == 3) return(c(0, 0,0,1,0,0))
            if (i == 4) return(c(1, 0,0,0,0,0))
            if (i == 5) return(c(0, 0,1,0,0,0))
            if (i == 6) return(c(0, 0,0,1,0,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))   
  
}