function y = code4(A,C,x0,M1,M2)

test = -1;

delta = 0.001;
N = 1000;
t = linspace(0,1,N+1);
h = 1/N;
h2 = h/2;

u = zeros(1,N+1);

x = zeros(1,N+1);
x(1) = x0;
lambda = zeros(1,N+1);

while(test < 0)
    
    oldu = u;
    oldx = x;
    oldlambda = lambda;
    
    for i = 1:N
        k1 = -0.5*x(i)^2 + C*u(i);
        k2 = -0.5*(x(i) + h2*k1)^2 + C*0.5*(u(i) + u(i+1));
        k3 = -0.5*(x(i) + h2*k2)^2 + C*0.5*(u(i) + u(i+1));
        k4 = -0.5*(x(i) + h*k3)^2 + C*u(i+1);
        x(i+1) = x(i) + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
    end
    
    for i = 1:N
        j = N + 2 - i;
        k1 = -A + lambda(j)*x(j);
        k2 = -A + (lambda(j) - h2*k1)*0.5*(x(j)+x(j-1));
        k3 = -A + (lambda(j) - h2*k2)*0.5*(x(j)+x(j-1));
        k4 = -A + (lambda(j) - h*k3)*x(j-1);
        lambda(j-1) = lambda(j) - (h/6)*(k1 + 2*k2 + 2*k3 + k4);
    end
    
    u1 = max(M1, min(C*lambda/2,M2));
    u = 0.5*(u1 + oldu);
        
    temp1 = delta*sum(abs(u)) - sum(abs(oldu - u));
    temp2 = delta*sum(abs(x)) - sum(abs(oldx - x));
    temp3 = delta*sum(abs(lambda)) - sum(abs(oldlambda - lambda));
    test = min(temp1, min(temp2, temp3));
end

y(1,:) = t;
y(2,:) = x;
y(3,:) = lambda;
y(4,:) = u;