function m = labten(p1,p2,p3,p4,p5,p6,p7)
global a b c x10 beta level T

a=p1;
b=p2;
c=p3;
x10=p4;
beta=p5;
level=p6;
T=p7;

solinit = bvpinit(linspace(0,T),@labteninit); 

sol = bvp4c(@labtenode,@labtenbc,solinit);

t=sol.x;
x1=sol.y(1,:);
x2=sol.y(2,:);
a1=sol.y(3,:);
a2=sol.y(4,:);
u=-(0.5/beta)*a2;
m(1,:)=t;
m(2,:)=x1;
m(3,:)=x2;
m(4,:)=u;

%--------------------------------------------------------------------------
function dydt=labtenode(t,y)
global a b c x10 beta level T

dydt=[-a*y(1) - b*y(2)
      -c*y(2) - (0.5/beta)*y(3)
      -2*(y(1) - level) + a*y(3)
      b*y(3) + c*y(4)];


%--------------------------------------------------------------------------
function res=labtenbc(ya,yb)
global a b c x10 beta level T

res=[ya(1)-x10
     ya(2)
     yb(3)    
     yb(4)];

%--------------------------------------------------------------------------
function v=labteninit(t)
global a b c x10 beta level T

v=[x10
   0.0
   0.0 
   0.0];
