function [w1,w2,u1,u2,yA,MH,k,x] = NonlineerAtis(a,b,alpha,beta,N,TOL,M)

% y'' = f(x,y,y') y(a)=alpha, y(b)=beta
% 1. y'' = f(x,y,y') y(a)=alpha y'(a)=t0 
% 2. z'' = f_y(x,y,y')z + f_y'(x,y,y')z' , z(a) = 0, z'(a)=1

% y = Y0, y' = Y1

f = @(x,Y0,Y1) ...;
fY0 = @(x,Y0,Y1) ...;
fY1 = @(x,Y0,Y1) ...;

h = (b-a)/N;

x = a:h:b;

t = (beta-alpha)/(b-a);

k = 1; 
while 1
    
    w1(1) = alpha;
    w2(1) = t; 
    
    for i = 1:N
        
        k11 = h*w2(i);
        k12 = h*f(x(i),w1(i),w2(i));
        
        k21 = h*(w2(i) + k12/2);
        k22 = h*f(x(i) + h/2, w1(i) + k11/2, w2(i) + k12/2);
        
        k31 = h*(w2(i) + k22/2);
        k32 = h*f(x(i) + h/2, w1(i) + k21/2, w2(i) + k22/2);
        
        k41 = h*(w2(i) + k32);
        k42 = h*f(x(i) + h, w1(i) + k31, w2(i) + k32);
        
        w1(i+1) = w1(i) + (1/6)*(k11 + 2*k21 + 2*k31 + k41);
        w2(i+1) = w2(i) + (1/6)*(k12 + 2*k22 + 2*k32 + k42);
        
    end
    
    plot(x,w1,'b-');
    hold on;
    
    if abs(w1(N+1) - beta) < TOL || k>M
        break;
    end
    
    u1(1) = 0;
    u2(1) = 1; 
    
    for i = 1:N
        
        m11 = h*u2(i);
        m12 = h*(fY0(x(i),w1(i),w2(i))*u1(i) + fY1(x(i),w1(i),w2(i))*u2(i));
        
        m21 = h*(u2(i) + m12/2);
        m22 = h*(fY0(x(i)+h/2,w1(i),w2(i))*(u1(i) + m11/2) + fY1(x(i)+h/2,w1(i),w2(i))*(u2(i) + m12/2));
        
        m31 = h*(u2(i) + m22/2);
        m32 = h*(fY0(x(i)+h/2,w1(i),w2(i))*(u1(i) + m21/2) + fY1(x(i)+h/2,w1(i),w2(i))*(u2(i) + m22/2));
        
        m41 = h*(u2(i) + m32);
        m42 = h*(fY0(x(i)+h,w1(i),w2(i))*(u1(i) + m31) + fY1(x(i)+h,w1(i),w2(i))*(u2(i) + m32));
        
        u1(i+1) = u1(i) + (1/6)*(m11 + 2*m21 + 2*m31 + m41);
        u2(i+1) = u2(i) + (1/6)*(m12 + 2*m22 + 2*m32 + m42);
        
    end
    
    t = t - (w1(N+1) - beta)/u1(N+1);
    
    k = k+1;
    
end

% Analitik çözüm

yA = ...;

% MutlakHata

MH = abs(yA - w1);