The problem 1 describes a uniformy loaded beam of variable stifness, simply supported at both end.

It may be defined as follow :

Solve the fourth order differential equation:

Subjected to the boundary conditions:

The exact solution of this problem is known to be:

 
N=1;// just one differential equation
m=4;//a fourth order  differential equation
M=sum(m);

x_low=1;x_up=2; // the x limits
zeta=[x_low,x_low,x_up,x_up]; //two constraints (on the value of u and its second derivative) on each bound.

//The external functions
//Theses functions are called by the solver with zu=[u(x);u'(x);u''(x);u'''(x)]

// - The function which computes the right hand side of the differential equation
function f=fsub(x,zu),f=(1-6*x^2*zu(4)-6*x*zu(3))/x^3,endfunction

// - The function which computes the derivative of fsub with respect to zu
function df=dfsub(x,zu),df=[0,0,-6/x^2,-6/x],endfunction

// - The function which computes the ith constraint for a given i
function g=gsub(i,zu),
  select i
  case 1 then  //x=zeta(1)=1
    g=zu(1) //u(1)=0
  case 2 then //x=zeta(2)=1
    g=zu(3) //u''(1)=0
  case 3 then //x=zeta(3)=2
    g=zu(1) //u(2)=0
  case 4 then  //x=zeta(4)=2
    g=zu(3) //u''(2)=0
  end
endfunction

// - The function which computes the derivative of gsub with respect to z
function dg=dgsub(i,z)
  select i
  case 1 then  //x=zeta(1)=1
    dg=[1,0,0,0]
  case 2 then //x=zeta(2)=1
    dg=[0,0,1,0]
  case 3 then //x=zeta(3)=2
     dg=[1,0,0,0]
  case 4 then  //x=zeta(4)=2
    dg=[0,0,1,0]
  end
endfunction

// - The function which computes the initial guess, unused here
function [zu,mpar]=guess(x),zu=0;mpar=0,endfunction 

 //define the function which computes the exact value of u for a given x ( for testing purposes)
function zu=trusol(x)
  zu=0*ones(4,1)
  zu(1) =  0.25*(10*log(2)-3)*(1-x) + 0.5 *( 1/x   + (3+x)*log(x) - x)
  zu(2) = -0.25*(10*log(2)-3)       + 0.5 *(-1/x^2 + (3+x)/x      + log(x) - 1)
  zu(3) = 0.5*( 2/x^3 + 1/x   - 3/x^2)
  zu(4) = 0.5*(-6/x^4 - 1/x/x + 6/x^3)
endfunction

fixpnt=[ ];//All boundary conditions are located at x_low and x_up

//    nonlin  collpnt n ntol ndimf  ndimi iprint iread iguess rstart nfxpnt
ipar=[0       0       1 2    2000   200   1      0     0      0      0     ]

ltol=[1,3];//set tolerance control on zu(1) and zu(3)
tol=[1.e-11,1.e-11];//set tolreance values for these two controls
xpoints=x_low:0.01:x_up;

zu=bvode(xpoints,N,m,x_low,x_up,zeta,ipar,ltol,tol,fixpnt,...
        fsub,dfsub,gsub,dgsub,guess)
//check the constraints
zu([1,3],[1 $]) //should be zero
plot(xpoints,zu(1,:)) // the evolution of the solution u
zu1=[];for x=xpoints,zu1=[zu1,trusol(x)]; end;  
norm(zu-zu1)
 

Same problem using bvodeS and an initial guess

 
function [z,lhS]=zstart(x)
  z=zeros(5,1);z(5)=1;
  lhS=[0;1];
endfunction
zu=bvode(xpoints,N,m,x_low,x_up,zeta,ltol=ltol,tol=tol,guess=zstart)
 

The problem 2 describes the small finite deformation of a thin shallow spherical cap of constant thickness subject to a quadratically varying axisymmetric external pressure distribution. Here phi is the meridian angle change of the deformed shell and psi is a stress function. For eps=mu=1e-3 two different solutions may found depending on the starting point

Subject to the boundary conditions

for x=0 and x=1

 
N=2;// two differential equations
m=[2 2];//each differential equation is of second  order
M=sum(m);

x_low=0;x_up=1; // the x limits
zeta=[x_low,x_low, x_up x_up]; //two  constraints on each bound.

//The external functions
//Theses functions are called by the solver with zu=[u1(x);u1'(x);u2(x);u2'(x)]

// - The function which computes the right hand side of the differential equation
function f=fsub2(x,zu,eps,dmu,eps4mu,gam,xt),
   f=[zu(1)/x^2-zu(2)/x+(zu(1)-zu(3)*(1-zu(1)/x)-gam*x*(1-x^2/2))/eps4mu //phi''
      zu(3)/x^2-zu(4)/x+zu(1)*(1-zu(1)/(2*x))/dmu];//psi''
endfunction

// - The function which computes the derivative of fsub with respect to zu
function df=dfsub2(x,zu,eps,dmu,eps4mu,gam,xt),
  df=[1/x^2+(1+zu(3)/x)/eps4mu, -1/x, -(1-zu(1)/x)/eps4mu, 0
      (1-zu(1)/x)/dmu             0    1/x^2              -1/x];
endfunction

// - The function which computes the ith constraint for a given i
function g=gsub2(i,zu),
  select i
  case 1 then  //x=zeta(1)=0
    g=zu(1) //u(0)=0
  case 2 then //x=zeta(2)=0
    g=-0.3*zu(3) //x*psi'-0.3*psi+0.7x=0
  case 3 then //x=zeta(3)=1
    g=zu(1) //u(1)=0
  case 4 then  //x=zeta(4)=1
    g=1*zu(4)-0.3*zu(3)+0.7*1 //x*psi'-0.3*psi+0.7x=0
  end
endfunction

// - The function which computes the derivative of gsub with respect to z
function dg=dgsub2(i,z)
  select i
  case 1 then  //x=zeta(1)=1
    dg=[1,0,0,0]
  case 2 then //x=zeta(2)=1
    dg=[0,0,-0.3,0]
  case 3 then //x=zeta(3)=2
     dg=[1,0,0,0]
  case 4 then  //x=zeta(4)=2
    dg=[0,0,-0.3,1]
  end
endfunction

gam=1.1
eps=1d-3
dmu=eps
eps4mu=eps^4/dmu
xt=sqrt(2*(gam-1)/gam)

fixpnt=[ ];//All boundary conditions are located at x_low and x_up
collpnt=4;
nsizef=4+3*M+(5+collpnt*N)*(collpnt*N+M)+(2*M-2)*2*M ;
nsizei=3 + collpnt*N+M;;
nmax=200;
//    nonlin  collpnt n  ntol ndimf        ndimi        iprint iread iguess rstart nfxpnt
ipar=[1       k       10 4    nmax*nsizef  nmax*nsizei   -1      0     0      0      0     ]

ltol=1:4;//set tolerance control on zu(1) zu(2 ) zu(3) and zu(4)
tol=[1.e-5,1.e-5,1.e-5,1.e-5];//set tolreance values for these four controls
xpoints=x_low:0.01:x_up;

zu=bvode(xpoints,N,m,x_low,x_up,zeta,ipar,ltol,tol,fixpnt,...
        fsub2,dfsub2,gsub2,dgsub2,guess2);
scf(1);clf();plot(xpoints,zu([1 3],:)) // the evolution of the solution phi and psi

//using an initial guess
// - The function which computes the initial guess, unused here
function [zu,dmval]=guess2(x,gam),
   cons=gam*x*(1-x^2/2)
   dcons=gam*(1-3*x^2/2)
   d2cons=-3*gam*x
   dmval=zeros(2,1)
   if x>xt then
     zu=[0 0 -cons -dcons]
     dmval(2)=-d2cons
   else
     zu=[2*x;2;-2*x+cons;-2*dcons]
     dmval(2)=d2cons
   end
endfunction 
ipar(9)=1;//iguess

zu2=bvode(xpoints,N,m,x_low,x_up,zeta,ipar,ltol,tol,fixpnt,...
        fsub2,dfsub2,gsub2,dgsub2,guess2);
scf(2);clf();plot(xpoints,zu2([1 3],:)) // the evolution of the solution phi and psi
 

An eigenvalue problem:

 
// y''(x)=-la*y(x)
// BV: y(0)=y'(0); y(1)=0
// Eigenfunctions and eigenvalues are y(x,n)=sin(s(n)*(1-x)), la(n)=s(n)^2,
// where s(n) are the zeros of f(s,n)=s+atan(s)-(n+1)*pi, n=0,1,2,...
// To get a third boundary condition, we choose y(0)=1
// (With y(x) also c*y(x) is a solution for each constant c.)
// We solve the following ode system:
// y''=-la*y
// la'=0
// BV: y(0)=y'(0), y(0)=1; y(1)=0
// z=[y(x) ; y'(x) ; la]

function rhs=fsub(x,z)
  rhs=[-z(3)*z(1);0]
endfunction

function g=gsub(i,z)
  g=[z(1)-z(2) z(1)-1 z(1)]
  g=g(i)
endfunction

// The following start function is good for the first 8 eigenfunctions.
function [z,lhs]=ystart(x,z,la0)
  z=[1;0;la0]
  lhs=[0;0]
endfunction

a=0;b=1;
m=[2;1];
n=2;
zeta=[a a b];
N=101;
x=linspace(a,b,N)';

// We have s(n)-(n+1/2)*pi -> 0 for n to infinity.
la0=input('n-th eigenvalue: n= ?');la0=(%pi/2+la0*%pi)^2;

z=bvodeS(x,m,n,a,b,fsub,gsub,zeta,ystart=list(ystart,la0));

clf()
plot(x,[z(1,:)' z(2,:)']) 
xtitle(['Startvalue =  '+string(la0);'Eigenvalue = '+string(z(3,1))],'x',' ')
legend(['y(x)';'y''(x)'])
 

A boundary value problem with more than one solution.

 
// DE: y''(x)=-exp(y(x))
// BV: y(0)=0; y(1)=0
// This boundary value problem has more than one solution.
// It is demonstrated how to find two of them with the help of
// some preinformation of the solutions y(x) to build the function ystart.
// z=[y(x);y'(x)]

a=0;b=1;m=2;n=1;
zeta=[a b];
N=101;
tol=1e-8*[1 1];
x=linspace(a,b,N);

function rhs=fsub(x,z),rhs=-exp(z(1));endfunction

function g=gsub(i,z)
  g=[z(1) z(1)]
  g=g(i)
endfunction

function [z,lhs]=ystart(x,z,M) 
  //z=[4*x*(1-x)*M ; 4*(1-2*x)*M]
  z=[M;0]
  //lhs=[-exp(4*x*(1-x)*M)]
  lhs=0
endfunction

for M=[1 4]
   if M==1
      z=bvodeS(x,m,n,a,b,fsub,gsub,zeta,ystart=list(ystart,M),tol=tol);
   else
      z1=bvodeS(x,m,n,a,b,fsub,gsub,zeta,ystart=list(ystart,M),tol=tol);
   end
end

// Integrating the ode yield e.g. the two solutions yex and yex1. 

function y=f(c),y=c.*(1-tanh(sqrt(c)/4).^2)-2;endfunction 
c=fsolve(2,f);

function y=yex(x,c)
  y=log(c/2*(1-tanh(sqrt(c)*(1/4-x/2)).^2))
endfunction 

function y=f1(c1), y=2*c1^2+tanh(1/4/c1)^2-1;endfunction
c1=fsolve(0.1,f1);

function y=yex1(x,c1)
  y=log((1-tanh((2*x-1)/4/c1).^2)/2/c1/c1)
endfunction 

disp(norm(z(1,:)-yex(x)),'norm(yex(x)-z(1,:))= ')
disp(norm(z1(1,:)-yex1(x)),'norm(yex1(x)-z1(1,:))= ')
clf();
subplot(2,1,1)
plot2d(x,z(1,:),style=[5])
xtitle('Two different solutions','x',' ') 
subplot(2,1,2)
plot2d(x,z1(1,:),style=[5])
xtitle(' ','x',' ')
 

A multi-point boundary value problem.

 
// DE y'''(x)=1
// z=[y(x);y'(x);y''(x)]
// BV: y(-1)=2 y(1)=2
// Side condition: y(0)=1 

a=-1;b=1;c=0;
// The side condition point c must be included in the array fixpnt.
n=1;
m=[3];

function rhs=fsub(x,z)
  rhs=1
endfunction

function g=gsub(i,z)
  g=[z(1)-2 z(1)-1 z(1)-2]
  g=g(i)
endfunction

N=10;
zeta=[a c b];
x=linspace(a,b,N);

z=bvodeS(x,m,n,a,b,fsub,gsub,zeta,fixpnt=c);
          
function y=yex(x)
y=x.^3/6+x.^2-x./6+1
endfunction

disp(norm(yex(x)-z(1,:)),'norm(yex(x)-z(1,:))= ')