function [wi, win, y, ti] = back_euler_1 ( t0, x0, tf, N, TOL, Nmax )
% Backward Euler con y'=-alpha(y-t^2)+2t
%   Puntode arranque, t0, x0; final, tf; punots N
 if nargin<6, Nmax=100; end, if nargin<5, TOL=0.00001; end, close all
 ti = [t0 zeros( 1,N )]; wi = zeros(1,N+1); win = zeros(1,N+1);
 wi(1) = x0'; win(1) = x0'; x01 = x0; y(1) = x0;
 h = ( tf-t0)/N;
 for i = 1:N
    w0 = x0;
    for j = 1:Nmax
       top = (w0 - x0) - h*beuler_1(t0+h,w0);
       bot = 1 - h*dbeuler_1(t0+h,w0);
       dw  = top/bot;
       w0  = w0-dw;
       if abs(dw)<TOL, break; end;
    end;
    x0  = w0;
    x01 = x01+h*beuler_1(t0,x01);
    t0 = t0 + h; ti(i+1) = t0;  
    wi(i+1) = x0'; win(i+1) = x01; y(i+1) = y(1)*exp(-10*t0)+t0^2; 
 end
 plot(ti,wi,ti,win,ti,y), xlabel('Tiempo'); ylabel('y')
 legend('Back. Euler', 'Euler', 'Exacto')
end

function f = beuler_1(t, p)
 f = -10*(p-t^2)+2*t;
end

function df = dbeuler_1(t, p)
 df = -10;
end