ex08ch3.m

来自「these codes are for solving OED with mat」· M 代码 · 共 22 行

function sol = ex08ch3
hold on;
for Z = [2 3 4]
    if Z == 2
        sol = bvpinit(linspace(0,Z,5),[0.5 -0.5]);
    else
        sol = bvpinit(sol,[0, Z]);
    end
    sol = bvp4c(@odes,@bcs,sol);
    plot(sol.x,sol.y(1,:),sol.x(end),sol.y(1,end),'ro')
    axis([0 4 0 1])
    drawnow
end
hold off
%======================================================
function dydz = odes(z,y)
alpha = 0.8;
dydz =  [ y(2); -2*z*y(2)/sqrt(1 - alpha*y(1)) ];

function res = bcs(ya,yb)
res = [ (ya(1) - 1) ; yb(1) ];