function hw5()

w = LinearFD()

s = NonlinearShooting()

function w = LinearFD( )

alpha = 1;
beta = 2;
xmin = 1;
xmax = 2;
h = 0.1;
x = xmin:h:xmax;

n = length(x);
w = zeros(n,1);
a = zeros(n,1);
b = zeros(n,1);
c = zeros(n,1);

w(1) = alpha;
w(n) = beta;

for i = 2:n-1;
    a(i) = 2+h^2*qx(x(i));
    b(i) = -1+h/2*px(x(i));
    c(i) = -1 - h/2*px(x(i));
    if (i == 2)
        d(i) = -h^2*rx(x(i))+(1+h/2*px(x(i)))*alpha;
    else if (i == n-1)
            d(i) = -h^2*rx(x(i))+(1-h/2*px(x(i)))*beta;
        else
            d(i) = -h^2*rx(x(i));
        end
    end   
end

for i = 3:n-1;
    a(i) = a(i)-c(i)/a(i-1)*b(i-1);
    d(i) = d(i)-c(i)/a(i-1)*d(i-1);
end

w(n-1) = d(n-1)/a(n-1);

for i = n-2:-1:2;
    w(i) = (d(i) - w(i+1)*b(i))/a(i);
end


function w1 = NonlinearShooting( )

xmin = 1;
xmax = 3;
alpha = 17;
beta = 43/3;
h = 0.1;
M = 1000;
k = 1;
TOL = 10^-5;

x = xmin:h:xmax;
n = length(x);
w1 = zeros(n,1);
w2 = zeros(n,1);

TK = (beta-alpha)/(xmax-xmin);

while (k <= M)
    w1(1) = alpha;
    w2(1) = TK;
    u1 = 0;
    u2 = 1;
    
    for i=2:n;
        k11 = h*w2(i-1);
        k21 = h*func(x(i-1),w1(i-1),w2(i-1));
        k12 = h*(w2(i-1)+1/2*k21);
        k22 = h*func(x(i-1)+h/2,w1(i-1)+1/2*k11,w2(i-1)+1/2*k21);
        k13 = h*(w2(i-1)+1/2*k22);
        k23 = h*func(x(i-1)+h/2,w1(i-1)+1/2*k12,w2(i-1)+1/2*k22);
        k14 = h*(w2(i-1)+k23);
        k24 = h*func(x(i),w1(i-1)+k13,w2(i-1)+k23);
        w1(i) = w1(i-1)+(k11+2*k12+2*k13+k14)/6;
        w2(i) = w2(i-1)+(k21+2*k22+2*k23+k24)/6;
        
        k11 = h*u2;
        k21 = h*(fy(x(i-1),w1(i-1),w2(i-1))*u1+fy1(x(i-1),w1(i-1),w2(i-1))*u2);
        k12 = h*(u2+1/2*k21);
        k22 = h*(fy(x(i-1)+h/2,w1(i-1),w2(i-1))*(u1+1/2*k11)+fy1(x(i-1)+h/2,w1(i-1),w2(i-1))*(u2+1/2*k21));
        k13 = h*(u2+1/2*k22);
        k23 = h*(fy(x(i-1)+h/2,w1(i-1),w2(i-1))*(u1+1/2*k12)+fy1(x(i-1)+h/2,w1(i-1),w2(i-1))*(u2+1/2*k22));
        k14 = h*(u2+k23);
        k24 = h*(fy(x(i),w1(i-1),w2(i-1))*(u1+k13)+fy1(x(i),w1(i-1),w2(i-1))*(u2+k23));
        u1 = u1+(k11+2*k12+2*k13+k14)/6;
        u2 = u2+(k21+2*k22+2*k23+k24)/6;
    end
    
    if (abs(w1(n)-beta) <= TOL)
        return;
    end
    TK = TK - (w1(n)-beta)/u1;
    k = k+1;
end

function y = func(x,y,y1)
y = 1/8*(32+2*x^3-y*y1);

function y = fy(x,y,y1)
y = -1/8*y1;

function y = fy1(x,y,y1)
y = -1/8*y;

function y = px(x)
y = -2/x;

function y = qx(x)
y = 2/(x^2);

function y = rx(x)
y = sin(log(x))/(x^2);



