Matlab - ayuda con matlab

Pues os comento,estoy haciendo el proyecto fin de carrera, estoy estudiando ingenieria química y tengo que hacer una simulacion de ozonizacion del fenol con matlab. Me han proporcionado una tesis que hace algo semejante para basarme,pero no consigo que funcione. La tesis reparte el trabajo en 4 archivos *.m. busco ayuda porque ya estoy desesperado. Incluso estaria dispuesto a dar una ayuda económica al que me ayude porque necesito acabarlo para septiembre como sea, ya que es lo unico que me queda.

EL primero expone el sistema de ecuaciones diferenciales que hay que resolver poniendolos en funcion de una variable y. es este:

function V=equation(y)
global Ha Ei R n;
Ha=2;Ei=10;R=100;
n=10;
h=1/n;
V=zeros(2*n,1);
y0=1;
V(1)=1/h^2*(-y0+2*y(1)-y(2))+Ha^2*y(1)*y(n+2);
for i=2:(n-1)
V(i)=1/h^2*(-y(i-1)+2*y(i)-y(i+1))+Ha^2*y(i)+(y(n+1+i));
end;
V(n)=(2/h^2)*(-y(n-1)+y(n)-(Ha^2-R)*h*y(n))+Ha^2*y(n);
V(n+1)=(2/h^2)*(y(n+1)-y(n+2))+Ha^2/(Ei-1)*y(n+1);
for i=(n+2):(2*n-1)
V(i)=1/h^2*(-y(i-1)+2*y(i)-y(i+1))+Ha^2/(Ei-1)*y(i-(n+1))*y(i);
end;
y2n1=1;
V(2*n)=1/h^2*(-y(2*n-1)+2*y(2*n)-y2n1)+Ha^2/(Ei-1)*y(n)*y(2*n);

el segundo archivo resuelve el sistema de ecuaciones igualando a cero, sacando los valores de y para averiguar dos variables que se usaran despues

global Ha Ei R n yi0 yi1;
Ha=2; Ei=10;R=100;n=10;
h=1/n;
X=[0:h:1];
yi0=ones(2*n,1);
options.iterations=1000;
options.MaxFunEvals=10000;
options.MaxIter=40000;
options.TolFun=1e-3;
options.TolX=1e-3;
y=fsolve('equation',yi0,options);
n=20;
h=1/n;
X1=[0:h:1];
yi1=interp1(X,y,X1);
options.TolFun=1e-5;
options.TolX=1e-5;
y=fsolve('equation',yi1,options);
A=[1;y(1:n)];
B=[x(n+1:2*n);1];
plot(X1,A,X1,B);
E=-1/(1-y(n))*(y(1)-1)/h;
D=-1/(1-y(n))*(y(n)-y(n-1))/h;
save tabla.mat E D

despues se escriben las ecuaciones correspondientes

function xdot=sistema(t,x,a)
global kl DA DB Cge m epsg kla nu k2;
a=kla/kl;
Qg=6e-5;
Vliq=8.5E-3;
Vgaz=epsg*Vliq
Ha=sqrt(k2*x(3)*DA/kl^2);
Ei=1+DB*x(3)/(nu*DA^(m*Cg));
R=k2*x(3)*(1-epsg)/kla;
Cg=(Cge-x(1))/log(Cge/x(1));
load(tabla.mat');
E=max(interp3(ha,ei,r,Etab,Ha,Ei,R),1);
D=min(interp3(ha,ei,r,Dtab,Ha,Ei,R),1);
xdot=zeros(3,1);
xdot(1)=(Qg+(Cge-x(1))-E*kla*(m*(Cge-x(1))/ln(Cge/x(1))-x(2))*Vliq)/(Vgaz*(1/ln(x(1)/Cge)+(Cge-x(1))/(x(1)*(ln(Cge/x(1))^2))))
xdot(2)=E*kla*(m*Cg-x(2))-((E-D)*kla*(m*Cg-x(2)))-k2*x(3)*x(2)*(1-a*DA/kl);
xdot(3)=-nu*(((E-D)*kla*(m*Cg-x(2)))+k2*x(3)*x(2)*(1-a*DA/kl));

y finalmente se resuelve:

global kl DA DB Cge m epsg kla nu k2;
kl=3.7e-4;
DA=1.7e-9;
DB=0.8e-9;
Cge=0.6;
m=0.3;
epsg=0.03;
kla=1.5e-2;
nu=1;
k2=10;
t0=0;
tf=4000;
CAgaz=1e-9;
CAliq=0;
CB=4.5;
x0=[CAgaz CAliq CB];
h=100;
[t,x]=ode23tb('sistema',[t0:tf/h:tf],x0,[]);
load(tabla.mat');
for i=1:h+1
Ha(i)=sqrt(k2*x(i,3)*DA/kl^2);
Ei(i)=1+DB*x(i,3)/(nu*DA*Cg(i)*m);
R(i)=k2*x(i,3)*(1-epsg)/kla;
Cg(i)=(Cge-x(i,1))/log(Cge/x(i,1));
E(i)=max(interp3(ha,ei,r,Etab,Ha1(i),Ei1(i),R(i)),1);
D(i)=min(interp3(ha,ei,r,Dtab,Ha1(i),Ei1(i),R(i)),1);
end;
figure(1);
subplot(2,2,1);
plot(t,Ha1);
title('Ha');
subplot(2,2,2);
plot(t,E);
title('E');
subplot(2,2,3);
plot(t,D);
title('D');
figure(2);
plot(t,x(:,1),'r');
figure(3);
plot(t,x(:,2),'r');
figure(4);
plot(t,x(:,3),'r');

eso es todo, y ya sabeis, si alguien esta interesado en ayudarme por unas pelillas aunque sea que se ponga en contacto conmigo en [email protected]

Un saludo.