function multisyn_ode()
display('Running multisyn_ode...');

%%%%%%%%%%%%%%%%%%%%%%%
%%%global parameters%%%
%%%%%%%%%%%%%%%%%%%%%%%

%%length of cable
L       =1000;

%%circumference of cable
l       =1;

%%diffusision coefficient
D       =0.1;

%%number of odes for U
stepx   =1;
x       =(0:stepx:L)';
M       =length(x);

%%utility zectors
o       =ones(size(x));
z       =zeros(size(x));

%%number and density of spines
rho     =1*o;%*(1:(2-1)/(M-1):2)';

%%%%%%%%%%%%%%%%%%%%%%
%%%local parameters%%%
%%%%%%%%%%%%%%%%%%%%%%

%%surface area of spine head without PSD
A       =1*o;%*(1:(2-1)/(M-1):2)';
% rgn     =rgn2idx([450 550],stepx);
% A(rgn)  =2*A(rgn);

%%surface area of PSD
a       =0.1*o;%*(1:(2-1)/(M-1):2)';
% rgn     =rgn2idx([450 550],stepx);
% Z(rgn)  =2*Z(rgn);

%%binding/unbinding rates
alpha   =0.0001*o;
beta    =0.0001*o;

%%scaffolding concentration
Z       =200*o;
% Z       =Z.*(1+1/2*sin(x/100));

%%spine neck hopping rate
omega   =0.000001*o;
% omega   =omega.*(1+1*sin(x/100));
% rgn         =rgn2idx([50 L],stepx);  
% omega(rgn)  =0*omega(rgn);
% omega       =0.001*(100*x+1);

%%PSD-ESM hopping rate
h       =0.001*o;
% rgn         =rgn2idx([400 600],stepx);  
% h(rgn)      =0.0000000001*h(rgn);

%%rate of endocytosis
k       =0.001*o;%0.001
if (0)
    rgn     =rgn2idx([90 110],stepx);
    k(rgn)  =10*k(rgn);
end

%%rate of exocytosis
sexo    =0.001*o;%*(1:(2-1)/(M-1):2)';%0.001
if (0)
    rgn         =rgn2idx([90 110],stepx);
    sexo(rgn)   =0.1*sexo(rgn);
end

%%rate of degradation
sdeg    =0.00001*o;%*(1:(2-1)/(M-1):2)';
if (0)
    rgn         =rgn2idx([90 110],stepx);
    sdeg(rgn)   =10*sdeg(rgn);
end

%%rate of production
delta   =0.001*o;%*(1:(2-1)/(M-1):2)';
if (0)
    rgn         =rgn2idx([90 110],stepx);
    delta(rgn)  =10*delta(rgn);
end

%%%%%%%%%%%%%%%%%%%%%%%%
%%%somatic parameters%%%
%%%%%%%%%%%%%%%%%%%%%%%%

%%rate of endocytosis
ksoma       =0.0001;

%%rate of exocytosis
ssoma       =0.0001;

%%rate of synthesis
deltasoma   =0.1;%0.1

%%rate of exit
kappa       =0.001;

%%flux from soma to dendrite
Jsoma       =deltasoma;

%colors
lgray  =[0.7 0.7 0.7];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%steady-state parameters%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

lambda  =sexo./(sexo+sdeg);
omegahat=omega.*k.*(1-lambda)./(omega+k.*(1-lambda));
Lambda2 =[stepx/2 stepx*ones(1,M-2) stepx/2].*(rho.*omegahat/D)';
r       =lambda.*delta./k./(1-lambda);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%calculation of steady-state%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I       =eye(M);
B       =G(x,x',L)*diag(Lambda2);
C       =I+B;
g       =G(x,0,L);
Cinvg   =C\g;
CinvB   =C\B;
Cinvo   =C\o;
if (Lambda2 == 0)
    if (omega == 0 && Jsoma ~= 0)
        error('No solution because omega == 0 but Jsoma ~= 0');
    end
    chi =0;
else
    chi =(Jsoma/D*(1-Lambda2*Cinvg)+Lambda2*r-Lambda2*CinvB*r)/...
        (Lambda2*Cinvo);
end
     
U0      =Jsoma/D*Cinvg+CinvB*r+chi*Cinvo;
R0      =(omega.*U0+lambda.*delta)./(omega+k.*(1-lambda));
P0      =(1+lambda.*k./h).*R0+lambda.*delta./h;
Q0      =alpha.*P0.*Z./(alpha.*P0+beta);
C0      =(k.*R0+delta)./(sexo+sdeg);
U0i     =zeros(size(U0));
R0i     =zeros(size(R0));
P0i     =zeros(size(P0));
Q0i     =zeros(size(Q0));
C0i     =zeros(size(C0));
Usoma0  =deltasoma/kappa;
Csoma0  =(ksoma+kappa)*deltasoma/kappa/ssoma;
Usoma0i =0;
Csoma0i =0;


%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%steady-state figures%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%

if(ishandle(1))
    clf(1);
end
figure(1);
p=plot(x,100*(U0+rho.*(A.*R0+C0+a.*(P0+Q0)))./...
    (U0(1)+rho(1)*(A(1)*R0(1)+C0(1)+a(1)*(P0(1)+Q0(1)))));
set(p,{'Color'},{'k'});
set(p,{'LineWidth'},{4});
set(p,{'LineStyle'},{'-'});
set(gca,'FontSize',16);
ylabel('receptors [% of soma]','FontWeight','bold');
xlabel('distance from soma [\mum]','FontWeight','bold');
% l=legend('total','GluR1/2','GluR2/3','Location','Best');
% set(l,'EdgeColor','white','Color','white');
% set(gca,'YLim',[20 100]);

if(ishandle(2))
    clf(2);
end
figure(2);
p=plot(x,[a.*(P0+Q0) a.*Q0]);
set(p,{'Color'},{'k'});
set(p,{'LineWidth'},{4;2});
set(p,{'LineStyle'},{'-'});
set(gca,'FontSize',16);
xlabel('distance from soma [\mum]','FontWeight','bold');
ylabel('PSD receptors [\mum^{-2} / spine]','FontWeight','bold');
l=legend('total','immobile','Location','Best');
set(l,'EdgeColor','white','Color','white','FontWeight','bold');
% if (delta1 == 0 && delta2 == 0)
%     set(gca,'YLim',[0 35]);
%     set(gca,'YTick',0:5:35);
% else
%     set(gca,'YLim',[0 60]);
%     set(gca,'YTick',0:10:60);
% end

if(ishandle(3))
    clf(3);
end
figure(3);
p=plot(x,100*Q0./(P0+Q0));
set(p,{'Color'},{'k'});
set(p,{'LineWidth'},{4});
set(p,{'LineStyle'},{'-'});
set(gca,'FontSize',16);
xlabel('distance from soma [\mum]','FontWeight','bold');
ylabel('immobile fraction [%]','FontWeight','bold');

if(ishandle(4))
    clf(4);
end
figure(4);
p=plot(x,A.*R0);
set(p,{'Color'},{'k'});
set(p,{'LineWidth'},{4});
set(p,{'LineStyle'},{'-'});
set(gca,'FontSize',16);
xlabel('distance from soma [\mum]','FontWeight','bold');
ylabel('ESM receptors [\mum^{-2} / spine]','FontWeight','bold');
% l=legend('total','GluR1/2','GluR2/3',...
%     'Location','Best');
% set(l,'EdgeColor','white','Color','white');
% if (delta1 == 0 && delta2 == 0)
%     set(gca,'YLim',[0 30]);
%     set(gca,'YTick',0:5:30);
% else
%     set(gca,'YLim',[0 80]);
%     set(gca,'YTick',0:10:80);
% end

if(ishandle(5))
    clf(5);
end
figure(5);
p=plot(x,C0);
set(p,{'Color'},{'k'});
set(p,{'LineWidth'},{4});
set(p,{'LineStyle'},{'-'});
set(gca,'FontSize',16);
xlabel('distance from soma [\mum]','FontWeight','bold');
ylabel('intracellular receptors [\mum^{-2} / spine]','FontWeight','bold');
% l=legend('total','GluR1/2','GluR2/3',...
%     'Location','Best');
% set(l,'EdgeColor','white','Color','white');
% if (delta1 == 0 && delta2 == 0)
%     set(gca,'YLim',[0 30]);
%     set(gca,'YTick',0:5:30);
% else
%     set(gca,'YLim',[0 80]);
%     set(gca,'YTick',0:10:80);
% end

if(ishandle(6))
    clf(6);
end
figure(6);
p=plot(x,U0);
set(p,{'Color'},{'k'});
set(p,{'LineWidth'},{4});
set(p,{'LineStyle'},{'-'});
set(gca,'FontSize',16);
xlabel('distance from soma [\mum]','FontWeight','bold');
ylabel('dendritic concentration [\mum^{-2}]','FontWeight','bold');
% l=legend('total','GluR1/2','GluR2/3','Location','Best');
% set(l,'EdgeColor','white','Color','white');
% if (delta1 == 0 && delta2 == 0)
%     set(gca,'YLim',[0 30]);
%     set(gca,'YTick',0:5:30);
% else
%     set(gca,'YLim',[0 80]);
%     set(gca,'YTick',0:10:80);
% end

if(ishandle(7))
    clf(7);
end
figure(7);
p=plot(x,[sexo.*C0 k.*A.*R0]);
set(p,{'Color'},{'k';'k'});
set(p,{'LineWidth'},{4;4});
set(p,{'LineStyle'},{'-';'--'});
set(gca,'FontSize',16);
l=legend('exocytosis','endocytosis','Location','Best');
set(l,'EdgeColor','white','Color','white','FontWeight','bold');
xlabel('distance from soma [\mum]','FontWeight','bold');
ylabel('rate [receptors \mum^{-2} s^{-1} / spine]','FontWeight','bold');

if(ishandle(8))
    clf(8);
end
figure(8);
p=plot(x/L,[U0 rho.*a.*(P0+Q0) rho.*a.*Q0 omega*10000]);
set(p,{'Color'},{'k';lgray;lgray;'r'});
set(p,{'LineWidth'},{4;4;2;2});
set(p,{'LineStyle'},{'-'});
set(gca,'FontSize',16);
% xlabel('distance from soma [\mum]','FontWeight','bold');
xlabel('distance from soma [mm]','FontWeight','bold');
ylabel('number or concentration [\mum^{-2}]','FontWeight','bold');
l=legend('concentration in dendrite','number in PSD','number bound in PSD',...
    'Location','Best');
set(l,'EdgeColor','white','Color','white','FontWeight','bold');
% set(gca,'XLim',[0 L]);
% set(gca,'XTick',0:50:L);
set(gca,'XLim',[0 1]);
set(gca,'XTick',0:0.1:1);
set(gca,'YLim',[0 220]);
set(gca,'YTick',0:20:220);
% if (delta1 == 0 && delta2 == 0)
%     set(gca,'YLim',[0 35]);
%     set(gca,'YTick',0:5:35);
% else
%     set(gca,'YLim',[0 60]);
%     set(gca,'YTick',0:10:60);
% end

% if(ishandle(8))
%     clf(8);
% end
% figure(8);
% p=plot(x,[U0 C0 rho.*a.*(P0+Q0) rho.*a.*Q0 40*o]);
% set(p,{'Color'},{'k';'k';lgray;lgray;lgray});
% set(p,{'LineWidth'},{4;2;4;2;4});
% set(p,{'LineStyle'},{'-';'-';'-';'-';'--'});
% set(gca,'FontSize',16);
% xlabel('distance from soma [\mum]','FontWeight','bold');
% ylabel('number or concentration [\mum^{-2}]','FontWeight','bold');
% % l=legend('concentration in dendrite','intracellular number',...
% %     'number in PSD','number bound in PSD','number in PSD before perturbation',...
% %     'Location','Best');
% % set(l,'EdgeColor','white','Color','white','FontWeight','bold');
% set(gca,'XLim',[0 200]);
% set(gca,'XTick',0:20:200);
% set(gca,'YLim',[0 80]);
% set(gca,'YTick',0:10:80);
% 
% if(ishandle(9))
%     clf(9);
% end
% figure(9);
% p=plot(x,100*a.*(P0+Q0)/40);
% set(p,{'Color'},{'k'});
% set(p,{'LineWidth'},{4});
% set(p,{'LineStyle'},{'-'});
% set(gca,'FontSize',16);
% xlabel('distance from soma [\mum]','FontWeight','bold');
% ylabel('PSD receptors [% of baseline]','FontWeight','bold');
% % l=legend('concentration in dendrite','intracellular number',...
% %     'number in PSD','number bound in PSD','number in PSD before perturbation',...
% %     'Location','Best');
% % set(l,'EdgeColor','white','Color','white','FontWeight','bold');
% set(gca,'XLim',[0 200]);
% set(gca,'XTick',0:20:200);
% set(gca,'YLim',[60 180]);
% set(gca,'YTick',60:10:180);

%%%scaling
% % a       =a.*(1+sin(x)/2);
% sexo    =sexo*2;
% k       =k*2;
% % delta   =delta*2;
% % Jsoma   =Jsoma*2;
% % alpha   =alpha.*(1+3/4*sin(x/10));
% % beta    =beta.*(1+3/4*cos(x/10));
% % Z   =Z.*(1+1/2*sin(x/10));
% 
% %initial conditions
% lambdaS  =sexo./(sexo+sdeg);
% omegahatS=omega.*k.*(1-lambdaS)./(omega+k.*(1-lambdaS));
% 
% Lambda2S =[stepx/2 stepx*ones(1,M-2) stepx/2].*(rho.*omegahatS/D)';
% rS       =lambdaS.*delta./k./(1-lambdaS);
% 
% BS       =G(x,x',L)*diag(Lambda2S);
% CS       =I+BS;
% CSinvg   =CS\g;
% CSinvBS  =CS\BS;
% CSinvo   =CS\o;
% if (Lambda2S == 0)
%     if (omega == 0 && Jsoma ~= 0)
%         error('No solution because omega == 0 but Jsoma ~= 0');
%     end
%     chi =0;
% else
%     chi =(Jsoma/D*(1-Lambda2S*CSinvg)+Lambda2S*rS-Lambda2S*CSinvBS*rS)/...
%         (Lambda2S*CSinvo);
% end
%      
% US0      =Jsoma/D*CSinvg+CSinvBS*rS+chi*CSinvo;
% RS0      =(omega.*US0+lambdaS.*delta)./(omega+k.*(1-lambdaS));
% PS0      =(1+lambdaS.*k./h).*RS0+lambdaS.*delta./h;
% QS0      =alpha.*PS0.*Z./(alpha.*PS0+beta);
% CS0      =(k.*RS0+delta)./(sexo+sdeg);
% UsomaS0  =deltasoma/kappa;
% CsomaS0  =(ksoma+kappa)*deltasoma/kappa/ssoma;
% 
% if(ishandle(9))
%     clf(9);
% end
% figure(9);
% p=plot(x/L,rho.*a.*(PS0+QS0));
% set(p,{'Color'},{'k'});
% set(p,{'LineWidth'},{4});
% set(p,{'LineStyle'},{'-'});
% set(gca,'FontSize',16);
% % xlabel('distance from soma [\mum]','FontWeight','bold');
% xlabel('distance from soma [mm]','FontWeight','bold');
% ylabel('PSD receptors [\mum^{-2}]','FontWeight','bold');
% % l=legend('baseline \sigma^{EXO}_0, k_0','Location','Best');
% l=legend('2\times \sigma^{EXO}_0, k_0','Location','Best');
% % l=legend('0.5\times \sigma^{EXO}_0, k_0','Location','Best');
% set(l,'EdgeColor','white','Color','white','FontWeight','bold');
% % set(gca,'XLim',[0 L]);
% % set(gca,'XTick',0:50:L);
% set(gca,'XLim',[0 1]);
% set(gca,'XTick',0:0.1:1);
% set(gca,'YLim',[0 50]);
% set(gca,'YTick',0:5:50);
% 
% if(ishandle(10))
%     clf(10);
% end
% figure(10);
% p=plot(x/L,100*((PS0+QS0)./(P0+Q0)-1));
% set(p,{'Color'},{'k'});
% set(p,{'LineWidth'},{4});
% set(p,{'LineStyle'},{'-'});
% set(gca,'FontSize',16);
% % xlabel('distance from soma [\mum]','FontWeight','bold');
% xlabel('distance from soma [mm]','FontWeight','bold');
% ylabel('PSD receptors [% change]','FontWeight','bold');
% % l=legend('concentration in dendrite','number in PSD','number bound in PSD',...
% %     'Location','Best');
% % set(l,'EdgeColor','white','Color','white','FontWeight','bold');
% % set(gca,'XLim',[0 L]);
% % set(gca,'XTick',0:50:L);
% set(gca,'XLim',[0 1]);
% set(gca,'XTick',0:0.1:1);
% set(gca,'YLim',[-30 50]);
% set(gca,'YTick',-30:10:50);
% % if (delta1 == 0 && delta2 == 0)
% %     set(gca,'YLim',[0 35]);
% %     set(gca,'YTick',0:5:35);
% % else
% %     set(gca,'YLim',[0 60]);
% %     set(gca,'YTick',0:10:60);
% % end

%%%%%%%%%%%%%%%%%%%%%
%%%simulation type%%%
%%%%%%%%%%%%%%%%%%%%%

kind=-1;
if (kind == 0)%test steady-state
    display('Testing multitraf on steady-state');
    tint    =[0 2*60*60];
elseif (kind == 1)%photoinactivation
    numhours    =24;
    display(['Photoinactivation for ' num2str(numhours) ' hours']);
    xplot       =floor([10 300])/stepx +1;
    recsoma     =0; %must be 0 (no) or 1 (yes)
    recdend     =0; %must be 0 (no) or 1 (yes)
    aprod       =1; %must be 0 (no) or 1 (yes)
    iprod       =0; %must be 0 (no) or 1 (yes)
    if (0)%no protein synthesis
        deltasoma   =0;
        delta       =0;
    end
elseif (kind == 2)%test steady-state
    display('Testing active multitrafi on steady-state');
    tint    =[0 2*60*60];
    recsoma     =0; %must be 0 (no) or 1 (yes)
    recdend     =0; %must be 0 (no) or 1 (yes)
    aprod       =1; %must be 0 (no) or 1 (yes)
    iprod       =0; %must be 0 (no) or 1 (yes)
elseif (kind == 3)%test steady-state
    display('Testing inactive multitrafi on steady-state');
    tint    =[0 2*60*60];
    recsoma     =1; %must be 0 (no) or 1 (yes)
    recdend     =1; %must be 0 (no) or 1 (yes)
    aprod       =0; %must be 0 (no) or 1 (yes)
    iprod       =1; %must be 0 (no) or 1 (yes)
elseif (kind == 4)%block exocytosis
    numminutes  =2*60;
    display(['Block exocytosis for ' num2str(numminutes) ' minutes']);
    xplot       =floor([100])/stepx +1;
    recsoma     =0; %must be 0 (no) or 1 (yes)
    recdend     =0; %must be 0 (no) or 1 (yes)
    aprod       =1; %must be 0 (no) or 1 (yes)
    iprod       =0; %must be 0 (no) or 1 (yes)
    ssoma       =0;
    sexo        =z;
elseif (kind == 5)%block endocytosis
    numminutes  =60;
    display(['Block endocytosis for ' num2str(numminutes) ' minutes']);
    xplot       =floor([100])/stepx +1;
    recsoma     =0; %must be 0 (no) or 1 (yes)
    recdend     =0; %must be 0 (no) or 1 (yes)
    aprod       =1; %must be 0 (no) or 1 (yes)
    iprod       =0; %must be 0 (no) or 1 (yes)
    ksoma       =0;
    k        =z;
end

%%%%%%%%%%%%%%%%%%%%%
%%%ode simulations%%%
%%%%%%%%%%%%%%%%%%%%%

Y   =[];
if (kind == 0)
    initY   =[U0;R0;P0;Q0;C0;Usoma0;Csoma0];
    [T,Y]   =ode45(@(t,y) multitraf(t,y,1,M,stepx,rho,A,a,Z,...
        D,alpha,beta,k,sexo,sdeg,delta,h,omega,ssoma,ksoma,kappa,deltasoma),tint,initY);
    
    U       =Y(:, 0*M+1: 1*M)';
    R       =Y(:, 1*M+1: 2*M)';
    P       =Y(:, 2*M+1: 3*M)';
    Q       =Y(:, 3*M+1: 4*M)';
    C       =Y(:, 4*M+1: 5*M)';
    Usoma   =Y(:, 5*M+1)';
    Csoma   =Y(:, 5*M+2)';
elseif (kind == 1)
    U       =z;
    R       =z;
    P       =z;
    Q       =z;
    C       =C0;
    Ui      =U0;
    Ri      =R0;
    Pi      =P0;
    Qi      =Q0;
    Ci      =z;

    Ux      =U(xplot);
    Rx      =R(xplot);
    Px      =P(xplot);
    Qx      =Q(xplot);
    Cx      =C(xplot);
    Uix     =Ui(xplot);
    Rix     =Ri(xplot);
    Pix     =Pi(xplot);
    Qix     =Qi(xplot);
    Cix     =Ci(xplot);
    
    Usoma   =0;
    Csoma   =Csoma0;
    
    Usomai  =Usoma0;
    Csomai  =0;

    T       =0;
    
    initY   =[U;R;P;Q;C;Ui;Ri;Pi;Qi;Ci;Usoma;Csoma;Usomai;Csomai];
    for i=1:numhours;
        display(['In hour ' num2str(i)]);
        [T0,Y]   =ode45(@(t,y) multitrafi(t,y,1,M,stepx,rho,A,a,Z,...
            D,alpha,beta,k,sexo,sdeg,delta,h,omega,ssoma,ksoma,kappa,deltasoma,...
            recsoma,recdend,aprod,iprod),[0 3600],initY);
        
        lenT0   =length(T0);
        U       =[U   Y(lenT0, 0*M+1: 1*M)'];
        R       =[R   Y(lenT0, 1*M+1: 2*M)'];
        P       =[P   Y(lenT0, 2*M+1: 3*M)'];
        Q       =[Q   Y(lenT0, 3*M+1: 4*M)'];
        C       =[C   Y(lenT0, 4*M+1: 5*M)'];
        Ui      =[Ui  Y(lenT0, 5*M+1: 6*M)'];
        Ri      =[Ri  Y(lenT0, 6*M+1: 7*M)'];
        Pi      =[Pi  Y(lenT0, 7*M+1: 8*M)'];
        Qi      =[Qi  Y(lenT0, 8*M+1: 9*M)'];
        Ci      =[Ci  Y(lenT0, 9*M+1:10*M)'];

        intT0   =2:lenT0;
        T       =[T;3600*(i-1)+T0(intT0)];
        Ux      =[Ux  Y(intT0, 0*M+xplot)'];
        Rx      =[Rx  Y(intT0, 1*M+xplot)'];
        Px      =[Px  Y(intT0, 2*M+xplot)'];
        Qx      =[Qx  Y(intT0, 3*M+xplot)'];
        Cx      =[Cx  Y(intT0, 4*M+xplot)'];
        Uix     =[Uix Y(intT0, 5*M+xplot)'];
        Rix     =[Rix Y(intT0, 6*M+xplot)'];
        Pix     =[Pix Y(intT0, 7*M+xplot)'];
        Qix     =[Qix Y(intT0, 8*M+xplot)'];
        Cix     =[Cix Y(intT0, 9*M+xplot)'];
        
        Usoma   =[Usoma  Y(intT0,10*M+1)'];
        Csoma   =[Csoma  Y(intT0,10*M+2)'];
        Usomai  =[Usomai Y(intT0,10*M+3)'];
        Csomai  =[Csomai Y(intT0,10*M+4)'];

        initY   =Y(lenT0,:)';        
    end
elseif (kind == 2)
    initY   =[U0;R0;P0;Q0;C0;z;z;z;z;z;Usoma0;Csoma0;0;0];
    [T,Y]   =ode45(@(t,y) multitrafi(t,y,1,M,stepx,rho,A,a,Z,...
        D,alpha,beta,k,sexo,sdeg,delta,h,omega,ssoma,ksoma,kappa,deltasoma,...
        recsoma,recdend,aprod,iprod),tint,initY);
    
    U       =Y(:, 0*M+1: 1*M)';
    R       =Y(:, 1*M+1: 2*M)';
    P       =Y(:, 2*M+1: 3*M)';
    Q       =Y(:, 3*M+1: 4*M)';
    C       =Y(:, 4*M+1: 5*M)';
    Usoma   =Y(:,10*M+1)';
    Csoma   =Y(:,10*M+2)';
elseif (kind == 3)
    initY   =[z;z;z;z;z;U0;R0;P0;Q0;C0;0;0;Usoma0;Csoma0];
    [T,Y]   =ode45(@(t,y) multitrafi(t,y,1,M,stepx,rho,A,a,Z,...
        D,alpha,beta,k,sexo,sdeg,delta,h,omega,ssoma,ksoma,kappa,deltasoma,...
        recsoma,recdend,aprod,iprod),tint,initY);
    
    U       =Y(:, 5*M+1: 6*M)';
    R       =Y(:, 6*M+1: 7*M)';
    P       =Y(:, 7*M+1: 8*M)';
    Q       =Y(:, 8*M+1: 9*M)';
    C       =Y(:, 9*M+1:10*M)';
    Usoma   =Y(:,10*M+3)';
    Csoma   =Y(:,10*M+4)';
elseif (kind == 4 || kind == 5)
    initY   =[U0;R0;P0;Q0;C0;Usoma0;Csoma0];
    [T,Y]   =ode45(@(t,y) multitraf(t,y,1,M,stepx,rho,A,a,Z,...
        D,alpha,beta,k,sexo,sdeg,delta,h,omega,ssoma,ksoma,kappa,deltasoma),...
        [0 numminutes*60],initY);
    
    U       =Y(:, 0*M+1: 1*M)';
    R       =Y(:, 1*M+1: 2*M)';
    P       =Y(:, 2*M+1: 3*M)';
    Q       =Y(:, 3*M+1: 4*M)';
    C       =Y(:, 4*M+1: 5*M)';
    Usoma   =Y(:, 5*M+1)';
    Csoma   =Y(:, 5*M+2)';
end

%%%%%%%%%%%%%%%%%%%%%%%%
%%%simulation figures%%%
%%%%%%%%%%%%%%%%%%%%%%%%

if (kind == 0 || kind == 2 || kind == 3)
    lenT    =length(T);
    Tplot   =[1 lenT];
    Aplot   =A*ones(size(Tplot));
    aplot   =a*ones(size(Tplot));

    if(ishandle(10))
        clf(10);
    end
    figure(10);
    p=plot(x,U(:,Tplot));
    set(p,{'Color'},{'k'});
    set(p,{'LineWidth'},{2});
    set(p,{'LineStyle'},{'-';'--'});
    xlabel('\bfdistance from soma [\mum]');
    ylabel('\bfconcentration in dendrite [\mum^{-2}]');

    if(ishandle(11))
        clf(11);
    end
    figure(11);
    p=plot(x,Aplot.*R(:,Tplot));
    set(p,{'Color'},{'k'});
    set(p,{'LineWidth'},{2});
    set(p,{'LineStyle'},{'-';'--'});
    xlabel('\bfdistance from soma [\mum]');
    ylabel('\bfreceptors in ESM');

    if(ishandle(12))
        clf(12);
    end
    figure(12);
    p=plot(x,C(:,Tplot));
    set(p,{'Color'},{'k'});
    set(p,{'LineWidth'},{2});
    set(p,{'LineStyle'},{'-';'--'});
    xlabel('\bfdistance from soma [\mum]');
    ylabel('\bfintracellular receptors');

    if(ishandle(13))
        clf(13);
    end
    figure(13);
    p=plot(x,aplot.*P(:,Tplot));
    set(p,{'Color'},{'k'});
    set(p,{'LineWidth'},{2});
    set(p,{'LineStyle'},{'-';'--'});
    xlabel('\bfdistance from soma [\mum]');
    ylabel('\bffree receptors in PSD');

    if(ishandle(14))
        clf(14);
    end
    figure(14);
    p=plot(x,aplot.*Q(:,Tplot));
    set(p,{'Color'},{'k'});
    set(p,{'LineWidth'},{2});
    set(p,{'LineStyle'},{'-';'--'});
    xlabel('\bfdistance from soma [\mum]');
    ylabel('\bfbound receptors in PSD');
    
elseif (kind == 1)
    titl    ='Photoinactivation';
    
    if (numhours < 3)
        numhours=numhours*60;
        T       =T/60;
        stime   ='min';
        steph   =10;
    elseif (numhours >= 3 && numhours < 24)
        T       =T/60/60;
        stime   ='hr';
        steph   =1;
    elseif (numhours >= 24 && numhours < 48)
        T       =T/60/60;
        stime   ='hr';
        steph   =2;
    elseif (numhours >= 48 && numhours < 96)
        T       =T/60/60;
        stime   ='hr';
        steph   =4;
    elseif (numhours >= 96)
        T       =T/60/60;
        stime   ='hr';
        steph   =8;
    end
    
    if(ishandle(10))
        clf(10);
    end
    figure(10);
    p=plot(T,Usoma);
    set(p,{'Color'},{'k'});
    set(p,{'LineWidth'},{4});
    set(p,{'LineStyle'},{'-'});
    set(gca,'FontSize',16);
%     l=legend('total','GluR1/2','GluR2/3','Location','Best');
%     set(l,'EdgeColor','white');
    title('soma');
    xlabel(['time [' stime ']'],'FontWeight','bold');
    ylabel('active surface receptors','FontWeight','bold');
    set(gca,'XLim',[0 numhours]);
    set(gca,'XTick',0:steph:numhours);
%     set(gca,'XLim',[0 3]);
%     set(gca,'XTick',0:0.5:3);

    if(ishandle(11))
        clf(11);
    end
    figure(11);
    p=plot(T,[ssoma.*Csoma;ksoma.*Usoma;kappa.*Usoma/l]);
    set(p,{'Color'},{'k';'k';'k'});
    set(p,{'LineWidth'},{4;4;2});
    set(p,{'LineStyle'},{'-';'--';'-'});
    set(gca,'FontSize',16);
    l=legend('exocytosis','endocytosis','flux','Location','Best');
    set(l,'EdgeColor','white','Color','white','FontWeight','bold');
    title('soma','FontWeight','bold');
    xlabel(['time [' stime ']'],'FontWeight','bold');
    ylabel('rate or flux [# of rec. s^{-1} or # of rec. s^{-1} \mum^{-1}]','FontWeight','bold');
    set(gca,'XLim',[0 numhours]);
    set(gca,'XTick',0:steph:numhours);
%     set(gca,'XLim',[0 3]);
%     set(gca,'XTick',0:0.5:3);

    if(ishandle(12))
        clf(12);
    end
    figure(12);
    p=plot(T,Csoma);
    set(p,{'Color'},{'k'});
    set(p,{'LineWidth'},{4});
    set(p,{'LineStyle'},{'-'});
    set(gca,'FontSize',16);
%     l=legend('total','GluR1/2','GluR2/3','Location','Best');
%     set(l,'EdgeColor','white');
    title('soma','FontWeight','bold');
    xlabel(['time [' stime ']'],'FontWeight','bold');
    ylabel('active intracellular receptors','FontWeight','bold');
    set(gca,'XLim',[0 numhours]);
    set(gca,'XTick',0:steph:numhours);
%     set(gca,'XLim',[0 3]);
%     set(gca,'XTick',0:0.5:3);

    if(ishandle(13))
        clf(13);
    end
    figure(13);
    p=plot(T,100*Usoma./Usoma0);
    set(p,{'Color'},{'k'});
    set(p,{'LineWidth'},{4});
    set(p,{'LineStyle'},{'-'});
    set(gca,'FontSize',16);
%     l=legend('total','GluR1/2','GluR2/3','Location','Best');
%     set(l,'EdgeColor','white');
    title('soma','FontWeight','bold');
    xlabel(['time [' stime ']'],'FontWeight','bold');
    ylabel('surface receptors [% of baseline]','FontWeight','bold');
    set(gca,'YLim',[0 100]);
    set(gca,'YTick',0:20:100);
    set(gca,'XLim',[0 numhours]);
    set(gca,'XTick',0:steph:numhours);
%     set(gca,'XLim',[0 3]);
%     set(gca,'XTick',0:0.5:3);

    for j=1:length(xplot)
        if(ishandle(4+10*j))
            clf(4+10*j);
        end
        figure(4+10*j);
        p=plot(T,a(xplot(j))*[Px(j,:)+Qx(j,:);Qx(j,:);Pix(j,:)+Qix(j,:);Qix(j,:)]);
        set(p,{'Color'},{'k';'k';lgray;lgray});
        set(p,{'LineWidth'},{4;2;4;2});
        set(p,{'LineStyle'},{'-'});
        set(gca,'FontSize',16);
        l=legend('active','immobile active','inactive','immobile inactive',...
            'Location','Best');
        set(l,'EdgeColor','white','Color','white','FontWeight','bold');
        title(['x = ' num2str((xplot(j)-1)*stepx) ' \mum'],'FontWeight','bold');
        xlabel(['time [' stime ']'],'FontWeight','bold');
        ylabel('PSD receptors [\mum^{-2} / spine]','FontWeight','bold');
        set(gca,'XLim',[0 numhours]);
        set(gca,'XTick',0:steph:numhours);

        if(ishandle(4+10*j+1))
            clf(4+10*j+1);
        end
        figure(4+10*j+1);
        p=plot(T,[sexo(xplot(j))*Cx(j,:);k(xplot(j))*A(xplot(j))*Rx(j,:)]);
        set(p,{'Color'},{'k';'k'});
        set(p,{'LineWidth'},{4;4});
        set(p,{'LineStyle'},{'-';'--'});
        set(gca,'FontSize',16);
        l=legend('exocytosis','endocytosis','Location','Best');
        set(l,'EdgeColor','white','Color','white','FontWeight','bold');
        title(['x = ' num2str((xplot(j)-1)*stepx) ' \mum'],'FontWeight','bold');
        xlabel(['time [' stime ']'],'FontWeight','bold');
        ylabel('rate [# of rec. s^{-1}]','FontWeight','bold');
        set(gca,'XLim',[0 numhours]);
        set(gca,'XTick',0:steph:numhours);

        if(ishandle(4+10*j+2))
            clf(4+10*j+2);
        end
        figure(4+10*j+2);
        p=plot(T,A(xplot(j))*[Rx(j,:);Rix(j,:)]);
        set(p,{'Color'},{'k';lgray});
        set(p,{'LineWidth'},{4});
        set(p,{'LineStyle'},{'-'});
        set(gca,'FontSize',16);
        l=legend('active','inactive','Location','Best');
        set(l,'EdgeColor','white','Color','white','FontWeight','bold');
        title(['x = ' num2str((xplot(j)-1)*stepx) ' \mum'],'FontWeight','bold');
        xlabel(['time [' stime ']'],'FontWeight','bold');
        ylabel('ESM receptors [\mum^{-2} / spine]','FontWeight','bold');
        set(gca,'XLim',[0 numhours]);
        set(gca,'XTick',0:steph:numhours);

        if(ishandle(4+10*j+3))
            clf(4+10*j+3);
        end
        figure(4+10*j+3);
        p=plot(T,Cx(j,:));
        set(p,{'Color'},{'k'});
        set(p,{'LineWidth'},{4});
        set(p,{'LineStyle'},{'-'});
        set(gca,'FontSize',16);
%         l=legend('total','GluR1/2','GluR2/3','Location','Best');
%         set(l,'EdgeColor','white','Color','white');
        title(['x = ' num2str((xplot(j)-1)*stepx) ' \mum'],'FontWeight','bold');
        xlabel(['time [' stime ']'],'FontWeight','bold');
        ylabel('active intracellular receptors [\mum^{-2} / spine]','FontWeight','bold');
        set(gca,'XLim',[0 numhours]);
        set(gca,'XTick',0:steph:numhours);

        if(ishandle(4+10*j+4))
            clf(4+10*j+4);
        end
        figure(4+10*j+4);
        p=plot(T,[Ux(j,:);Uix(j,:)]);
        set(p,{'Color'},{'k';lgray});
        set(p,{'LineWidth'},{4});
        set(p,{'LineStyle'},{'-'});
        set(gca,'FontSize',16);
        l=legend('active','inactive','Location','Best');
        set(l,'EdgeColor','white','Color','white','FontWeight','bold');
        title(['x = ' num2str((xplot(j)-1)*stepx) ' \mum'],'FontWeight','bold');
        xlabel(['time [' stime ']'],'FontWeight','bold');
        ylabel('dendritic concentration [\mum^{-2}]','FontWeight','bold');
        set(gca,'XLim',[0 numhours]);
        set(gca,'XTick',0:steph:numhours);

        if(ishandle(4+10*j+5))
            clf(4+10*j+5);
        end
        figure(4+10*j+5);
        p=plot(T,100*(Px(j,:)+Qx(j,:))/(P0(xplot(j))+Q0(xplot(j))));
        set(p,{'Color'},{'k'});
        set(p,{'LineWidth'},{4});
        set(p,{'LineStyle'},{'-'});
        set(gca,'FontSize',16);
%         l=legend('total','GluR1/2','GluR2/3','Location','Best');
%         set(l,'EdgeColor','white','Color','white');
        title(['x = ' num2str((xplot(j)-1)*stepx) ' \mum'],'FontWeight','bold');
        xlabel(['time [' stime ']'],'FontWeight','bold');
        ylabel('PSD receptors [% of baseline]','FontWeight','bold');
        set(gca,'XLim',[0 numhours]);
        set(gca,'XTick',0:steph:numhours);

        if(ishandle(4+10*j+6))
            clf(4+10*j+6);
        end
        figure(4+10*j+6);
        p=plot(T,100*Rx(j,:)/R0(xplot(j)));
        set(p,{'Color'},{'k'});
        set(p,{'LineWidth'},{4});
        set(p,{'LineStyle'},{'-'});
        set(gca,'FontSize',16);
%         l=legend('total','GluR1/2','GluR2/3','Location','Best');
%         set(l,'EdgeColor','white','Color','white');
        title(['x = ' num2str((xplot(j)-1)*stepx) ' \mum'],'FontWeight','bold');
        xlabel(['time [' stime ']'],'FontWeight','bold');
        ylabel('ESM receptors [% of baseline]','FontWeight','bold');
        set(gca,'XLim',[0 numhours]);
        set(gca,'XTick',0:steph:numhours);

        if(ishandle(4+10*j+7))
            clf(4+10*j+7);
        end
        figure(4+10*j+7);
        p=plot(T,100*(a(xplot(j))*(Px(j,:)+Qx(j,:))+A(xplot(j))*Rx(j,:))/...
            (a(xplot(j))*(P0(xplot(j))+Q0(xplot(j)))+A(xplot(j))*R0(xplot(j))));
        set(p,{'Color'},{'k'});
        set(p,{'LineWidth'},{4});
        set(p,{'LineStyle'},{'-'});
        set(gca,'FontSize',16);
%         l=legend('total','GluR1/2','GluR2/3','Location','Best');
%         set(l,'EdgeColor','white','Color','white');
        title(['x = ' num2str((xplot(j)-1)*stepx) ' \mum'],'FontWeight','bold');
        xlabel(['time [' stime ']'],'FontWeight','bold');
        ylabel('receptors in spines [% of baseline]','FontWeight','bold');
        set(gca,'XLim',[0 numhours]);
        set(gca,'XTick',0:steph:numhours);

        if(ishandle(4+10*j+8))
            clf(4+10*j+8);
        end
        figure(4+10*j+8);
        p=plot(T,100*Ux(j,:)/U0(xplot(j)));
        set(p,{'Color'},{'k'});
        set(p,{'LineWidth'},{4});
        set(p,{'LineStyle'},{'-'});
        set(gca,'FontSize',16);
%         l=legend('total','GluR1/2','GluR2/3','Location','Best');
%         set(l,'EdgeColor','white','Color','white');
        title(['x = ' num2str((xplot(j)-1)*stepx) ' \mum'],'FontWeight','bold');
        xlabel(['time [' stime ']'],'FontWeight','bold');
        ylabel('dendritic receptors [% of baseline]','FontWeight','bold');
        set(gca,'XLim',[0 numhours]);
        set(gca,'XTick',0:steph:numhours);

        if(ishandle(4+10*j+9))
            clf(4+10*j+9);
        end
        figure(4+10*j+9);
        p=plot(T,100*(rho(xplot(j))*(a(xplot(j))*(Px(j,:)+Qx(j,:))...
            +A(xplot(j))*Rx(j,:))+Ux(j,:))/...
            (rho(xplot(j))*(a(xplot(j))*(P0(xplot(j))+Q0(xplot(j)))...
            +A(xplot(j))*R0(xplot(j)))+U0(xplot(j))));
        set(p,{'Color'},{'k'});
        set(p,{'LineWidth'},{4});
        set(p,{'LineStyle'},{'-'});
        set(gca,'FontSize',16);
%         l=legend('total','GluR1/2','GluR2/3','Location','Best');
%         set(l,'EdgeColor','white','Color','white');
        title(['x = ' num2str((xplot(j)-1)*stepx) ' \mum'],'FontWeight','bold');
        xlabel(['time [' stime ']'],'FontWeight','bold');
        ylabel('surface receptors [% of baseline]','FontWeight','bold');
        set(gca,'XLim',[0 numhours]);
        set(gca,'XTick',0:steph:numhours);
    end
elseif (kind == 4 || kind == 5)
    T       =T/60;
    stime   ='min';
    steph   =10;
    
    for j=1:length(xplot)
        if(ishandle(5+j*5))
            clf(5+j*5);
        end
        figure(5+j*5);
        p=plot(T,a(xplot(j))*[P(xplot(j),:)+Q(xplot(j),:);Q(xplot(j),:)]);
        set(p,{'Color'},{'k'});
        set(p,{'LineWidth'},{4;2});
        set(p,{'LineStyle'},{'-'});
        set(gca,'FontSize',16,'FontWeight','bold');
        title(['x = ' num2str((xplot(j)-1)*stepx) ' \mum']);
        xlabel(['time [' stime ']']);
        ylabel('PSD receptors [\mum^{-2} / spine]');
        set(gca,'XLim',[0 numminutes]);
        set(gca,'XTick',0:steph:numminutes);

        if(ishandle(5+j*5+1))
            clf(5+j*5+1);
        end
        figure(5+j*5+1);
        p=plot(T,A(xplot(j))*R(xplot(j),:));
        set(p,{'Color'},{'k'});
        set(p,{'LineWidth'},{4});
        set(p,{'LineStyle'},{'-'});
        set(gca,'FontSize',16,'FontWeight','bold');
        title(['x = ' num2str((xplot(j)-1)*stepx) ' \mum']);
        xlabel(['time [' stime ']']);
        ylabel('ESM receptors [\mum^{-2} / spine]');
        set(gca,'XLim',[0 numminutes]);
        set(gca,'XTick',0:steph:numminutes);

        if(ishandle(5+j*5+2))
            clf(5+j*5+2);
        end
        figure(5+j*5+2);
        p=plot(T,C(xplot(j),:));
        set(p,{'Color'},{'k'});
        set(p,{'LineWidth'},{4});
        set(p,{'LineStyle'},{'-'});
        set(gca,'FontSize',16,'FontWeight','bold');
        title(['x = ' num2str((xplot(j)-1)*stepx) ' \mum']);
        xlabel(['time [' stime ']']);
        ylabel('intracellular receptors [\mum^{-2} / spine]');
        set(gca,'XLim',[0 numminutes]);
        set(gca,'XTick',0:steph:numminutes);

        if(ishandle(5+j*5+3))
            clf(5+j*5+3);
        end
        figure(5+j*5+3);
        p=plot(T,U(xplot(j),:));
        set(p,{'Color'},{'k'});
        set(p,{'LineWidth'},{4});
        set(p,{'LineStyle'},{'-'});
        set(gca,'FontSize',16,'FontWeight','bold');
        title(['x = ' num2str((xplot(j)-1)*stepx) ' \mum']);
        xlabel(['time [' stime ']']);
        ylabel('dendritic receptors [\mum^{-2}]');
        set(gca,'XLim',[0 numminutes]);
        set(gca,'XTick',0:steph:numminutes);

        if(ishandle(5+j*5+4))
            clf(5+j*5+4);
        end
        figure(5+j*5+4);
        p=plot(T,(P(xplot(j),:)+Q(xplot(j),:))/(P(xplot(j),1)+Q(xplot(j),1)));
        set(p,{'Color'},{'k'});
        set(p,{'LineWidth'},{4});
        set(p,{'LineStyle'},{'-'});
        set(gca,'FontSize',16,'FontWeight','bold');
        title(['x = ' num2str((xplot(j)-1)*stepx) ' \mum']);
        xlabel(['time [' stime ']']);
        ylabel('PSD receptors [\mum^{-2} / spine]');
        set(gca,'XLim',[0 numminutes]);
        set(gca,'XTick',0:steph:numminutes);
    end
end
clear all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%auxilary function: rgn2idx%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function idx=rgn2idx(rgn,stepx)
idx =floor(rgn(1)/stepx)+1:ceil(rgn(2)/stepx)+1;

%%%%%%%%%%%%%%%%%%%%%%
%%%Green's function%%%
%%%%%%%%%%%%%%%%%%%%%%

function z=G(r,c,L)
or  =ones(size(c));
oc  =ones(size(r));
P   =(r*or+oc*c)/L;
M   =(r*or-oc*c)/L;
z   =L/12*(h(P)+h(M));

function z=h(x)
z   =3*x.^2-6*abs(x)+2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%ode function for trafficking%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dy = multitraf(t,y,l,M,stepx,rho,A,a,Z,...
        D,alpha,beta,k,sexo,sdeg,delta,h,omega,ssoma,ksoma,kappa,deltasoma)
    
dy  =zeros(5*M+2,1);
U       = 0*M+1: 1*M;
R       = 1*M+1: 2*M;
P       = 2*M+1: 3*M;
Q       = 3*M+1: 4*M;
C       = 4*M+1: 5*M;
Usoma   = 5*M+1;
Csoma   = 5*M+2;

dy(U(1))    = 2/stepx*(kappa*y(Usoma)/l+D*(y(U(2))-y(U(1)))/stepx)...
             -rho(1)*omega(1)*(y(U(1))-y(R(1)));
dy(U(2:M-1))= D/(stepx^2)*(y(U(3:M))-2*y(U(2:M-1))+y(U(1:M-2)))...
             -rho(2:M-1).*omega(2:M-1).*(y(U(2:M-1))-y(R(2:M-1)));
dy(U(M))    =-2*D/(stepx^2)*(y(U(M))-y(U(M-1)))...
             -rho(M)*omega(M)*(y(U(M))-y(R(M)));
dy(R)       =(omega.*(y(U)-y(R))-k.*y(R)-h.*(y(R)-y(P)))./A;
dy(P)       =(h.*(y(R)-y(P))+sexo.*y(C))./a-alpha.*(Z-y(Q)).*y(P)+beta.*y(Q);
dy(Q)       = alpha.*(Z-y(Q)).*y(P)-beta.*y(Q);
dy(C)       =-(sexo+sdeg).*y(C)+k.*y(R)+delta;

dy(Usoma)   = ssoma*y(Csoma)-(ksoma+kappa)*y(Usoma);
dy(Csoma)   =-ssoma*y(Csoma)+ksoma*y(Usoma)+deltasoma;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%ode function for trafficking with inactivation%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dy = multitrafi(t,y,l,M,stepx,rho,A,a,Z,...
        D,alpha,beta,k,sexo,sdeg,delta,h,omega,ssoma,ksoma,kappa,deltasoma,...
        recsoma,recdend,aprod,iprod)
    
dy      =zeros(10*M+4,1);
U       = 0*M+1: 1*M;
R       = 1*M+1: 2*M;
P       = 2*M+1: 3*M;
Q       = 3*M+1: 4*M;
C       = 4*M+1: 5*M;
Ui      = 5*M+1: 6*M;
Ri      = 6*M+1: 7*M;
Pi      = 7*M+1: 8*M;
Qi      = 8*M+1: 9*M;
Ci      = 9*M+1:10*M;
Usoma   =10*M+1;
Csoma   =10*M+2;
Usomai  =10*M+3;
Csomai  =10*M+4;

dy(U(1))    = 2/stepx*(kappa*y(Usoma)/l+D*(y(U(2))-y(U(1)))/stepx)...
             -rho(1)*omega(1)*(y(U(1))-y(R(1)));
dy(U(2:M-1))= D/(stepx^2)*(y(U(3:M))-2*y(U(2:M-1))+y(U(1:M-2)))...
             -rho(2:M-1).*omega(2:M-1).*(y(U(2:M-1))-y(R(2:M-1)));
dy(U(M))    =-2*D/(stepx^2)*(y(U(M))-y(U(M-1)))...
             -rho(M)*omega(M)*(y(U(M))-y(R(M)));
dy(R)       =(omega.*(y(U)-y(R))-k.*y(R)-h.*(y(R)-y(P)))./A;
dy(P)       =(h.*(y(R)-y(P))+sexo.*y(C))./a-alpha.*(Z-y(Q)-y(Qi)).*y(P)+beta.*y(Q);
dy(Q)       = alpha.*(Z-y(Q)-y(Qi)).*y(P)-beta.*y(Q);
dy(C)       =-(sexo+sdeg).*y(C)+k.*y(R)+aprod*delta;

dy(Ui(1))    = 2/stepx*(kappa*y(Usomai)/l+D*(y(Ui(2))-y(Ui(1)))/stepx)...
             -rho(1)*omega(1)*(y(Ui(1))-y(Ri(1)));
dy(Ui(2:M-1))= D/(stepx^2)*(y(Ui(3:M))-2*y(Ui(2:M-1))+y(Ui(1:M-2)))...
             -rho(2:M-1).*omega(2:M-1).*(y(Ui(2:M-1))-y(Ri(2:M-1)));
dy(Ui(M))    =-2*D/(stepx^2)*(y(Ui(M))-y(Ui(M-1)))...
             -rho(M)*omega(M)*(y(Ui(M))-y(Ri(M)));
dy(Ri)       =(omega.*(y(Ui)-y(Ri))-k.*y(Ri)-h.*(y(Ri)-y(Pi)))./A;
dy(Pi)       =(h.*(y(Ri)-y(Pi))+recdend*sexo.*y(Ci))./a...
              -alpha.*(Z-y(Q)-y(Qi)).*y(Pi)+beta.*y(Qi);
dy(Qi)       = alpha.*(Z-y(Q)-y(Qi)).*y(Pi)-beta.*y(Qi);
dy(Ci)       =-(recdend*sexo+sdeg).*y(Ci)+k.*y(Ri)+iprod*delta;

dy(Usoma)   = ssoma*y(Csoma)-(ksoma+kappa)*y(Usoma);
dy(Csoma)   =-ssoma*y(Csoma)+ksoma*y(Usoma)+aprod*deltasoma;

dy(Usomai)   = recsoma*ssoma*y(Csomai)-(ksoma+kappa)*y(Usomai);
dy(Csomai)   =-recsoma*ssoma*y(Csomai)+ksoma*y(Usomai)+iprod*deltasoma;