Re: Problem with Lsode..help

[genehacker]

Hi Bryan, Here is the script

**********************

clf;
clc;
clear all;
1;

pre_OB = 1000;
pre_OC = 1;
f1 = 1.8; IC2 = 5;
p = 2;
k1min = 0; k1max =  0.9*p;
k2min = 0; k2max =  0.5*p;
k23min = 0; k23max =  0.03*p;
k3 = 1.0*p;               
f2 = 1.8; kI = 1; vmax = 0.01; 


global x4; 

x5 = @(x1)   k1max- ((k1max- k1min)/(1 + (f1 * x1/IC2)));
x6 = @(x1)   k2max- ((k2max- k2min)/(1 + (0.01*f1 * x1/IC2)));
x7 = @(x1)   k23min- ((k23min- k23max)/(1 + (0.01*f1 * x1/IC2)));
t0 = 8;

x8 = @(t,x2) ((t+0.001) .* (t<t0)) + (t0 .*(t==t0)) +
(t0*exp(0.1*t0)*exp(-0.1*t) .*(t > t0));
xdot = @(x,t) [ ( (x8(t,x(2)) * x4/pre_OC) - x5(x(1)) * x(1));
                        (x6(x(1)) * pre_OB - x7(x(1)) * x(2));
                        (x7(x(1)) * x(2) - k3 * x(3)) ];

x = [0,0,0; 0,0,0];
len = 100;
time = linspace(0,len,len+1);

matrix = zeros(1,10);
matrix(1,10) = bmd = 100; % for bone initial state
a = 0.06;    % OB weight
b = 0.1;     % OC weight
for i=1:(length(time)-1)
        seedtime =  time(i:i+1);
        [x,istate,msg] = lsode(xdot,[x(2,1),x(2,2),x(2,3)],seedtime);
        bmd = bmd + a*x(2,3) - b*x(2,1);
        x4 = (vmax * pre_OC) / ( x8(time(i+1),x(2,2)) * (pre_OC + (
x8(time(i+1),x(2,2)) * (1 + (10 * x(2,2) / kI))  )) );
        matrix(i+1,:) = [time(i+1) x8(time(i+1),x(2,2)) x(2,1) x(2,2) x(2,3)
x5(x(2,1)) x6(x(2,1)) x7(x(2,1)) x4 bmd];
endfor

timerescale = 1;
plot(time/timerescale,10*matrix(:,3),"r"); hold on ;
plot(time/timerescale,10*matrix(:,5),"m"); hold on;
plot(time/timerescale,matrix(:,10),"b"); hold off;
xlabel("DAYS");
ylabel('Concentration');
axis([0 len]);
legend('OC','OB','bmd');
print -djpg foo.jpg;

*************************

-- 
View this message in context: 
http://www.nabble.com/Problem-with-Lsode..help-tp20193330p20195750.html
Sent from the Octave - General mailing list archive at Nabble.com.