ex934.m

ch9例题05_ode23s口服药在人体中的吸收 matlab 程序

function ex934
clear;
e=0.5e-6; i=1; m(1)=0; m(2)=10; yn=35/3; h=0.02; N=1/h;
while 1
    x(1)=2; y(:,1)=[8;m(i)];
    for n=1:N
        k1=f_shoot(y(:,n),x(n));
        k2=f_shoot(y(:,n)+h*k1/2,x(n)+h/2);
        k3=f_shoot(y(:,n)+h*k2/2,x(n)+h/2);
        k4=f_shoot(y(:,n)+h*k3,x(n)+h);
        y(:,n+1)=y(:,n)+h*(k1+2*k2+2*k3+k4)/6;
        x(n+1)=x(n)+h;
    end
    y1_end(i)=y(1,n+1);
    if (abs(y1_end(i)-yn)<e) [x',y'], break, end
    if i>=2 
        m(i+1)=m(i)-(m(i)-m(i-1))*(y1_end(i)-yn)/(y1_end(i)-y1_end(i-1));
    end
    i=i+1;
end
m'

plot(x,y(1,:),'-r',x,y(2,:),'.-b');
xlabel('x'),ylabel('y(x) and y{\prime}(x)')
text(2.1,7,['Guess for y{\prime}(0)=',num2str(m(i))])
text(2.6,7,['y(3)=',num2str(y1_end(i))])
gtext('y(x)');gtext('y{\prime}(x)');

%----------------------------------------------------------------------
function f=f_shoot(y,x)
f=[y(2);-y(1)*y(2)/4+x^3/2+4];