%% Laboration 2 - Vecka 2
%% 2 Begynnelsevärdesproblem 
% Uppg 1
hold off
uprim=@(t,u) -u+sin(t)+cos(t); % högerled. dvs. u'(t)
f=@(t,u) sin(t)+u.*exp(-t);
a=0;b=4;n=20;
t=linspace(a,b,n); u=linspace(-2,2,n);
[T,U]=meshgrid(t,u);
DT=ones(size(T)), DU=uprim(T,U)
quiver(T,U,DT,DU,0.9), hold on, grid on

plot(t, f(t, 0), 'g--',t, f(t, -1), 'g--',t, f(t, 1), 'r-',t, f(t, 2), 'g--'  )

%% Uppg 2
hold off
uprim=@(t,u) cos(3.*t)-sin(5*t).*u;
a=0;b=5;N=20;
t=linspace(a,b,N); % För N intervall krävs N+1 punkter.
u=linspace(0, 3, N);
[T,U]=meshgrid(t,u);
DT=ones(size(T)); DU=uprim(T,U);
quiver(T,U,DT,DU, 0.8), hold on

for N=[10 30 50 100]
    h=(b-a)/N
    t=linspace(a,b,N+1);u=zeros(size(t))
    ua=2;
    u(1)=ua; % at u(t), t=0
    for n=1:N
        u(n+1)=u(n)+h*uprim(t(n), u(n))
    end
    plot(t, u)
end

%% 2.3 Konvergens
% Uppg 3
uprim=@(t,u) cos(3*t)-sin(5*t).*u
ua=2; %begynnelsevärde i u(0)
a=0; b=5; n=40;
t=linspace(a,b,n);u=linspace(0,3,n);
[T,U]=meshgrid(t,u);
DT=ones(size(T));DU=uprim(T,U);
quiver(T,U,DT,DU,0.9), hold on

h=(b-a)/n
u=zeros(size(t))
u(1)=ua; % at u(t), t=0
for n=1:n-1
    u(n+1)=u(n)+h*uprim(t(n), u(n))
end
plot(t, u), hold on

[t,u]=ode45(uprim,[a b], ua)
plot(t,u), grid on

legend('','Eulers framåtmetod', 'ode45')

%% 2.5 System av första ordningens differentialekvationer 
f1=@(u1,u2)2*u1-u1.*u2;
f2=@(u1,u2)u2+0.4*u1.*u2-u2.^2;
f=@(t,u)[f1(u(1),u(2)); f2(u(1),u(2))];
T=8; tspan=linspace(0,T,200);
u0=[0.1;0.1];
[t,U]=ode45(f,tspan,u0);
% plot(t,U(:,1), 'b', t, U(:,2), 'g'), hold on, grid on

% Fasporträtt
u1=linspace(-1,4,25); u2=linspace(-1,3,25);
[U1,U2]=meshgrid(u1,u2);
quiver(U1,U2,f1(U1,U2),f2(U1,U2), 1.5), hold on, grid on
axis([-1 5 -1 3]) %axis bounds
plot(U(:,1), U(:,2), 'r')
xlabel('u_1');ylabel('u_2');

plot(0,0,'r.', 0,1,'r.', 2.5 ,2,'r.')

%% Uppg 4
% Funktioner u'1 och u'2 i ekvationssystem
f1=@(u1,u2) cos(1+u1-u2.^2+0.3)
f2=@(u1,u2) sin(0.3+u1+u2-0.4*u1);
f=@(t,u)    [f1(u(1),u(2)); f2(u(1),u(2))];
% Begynnelsevärden u1(0) och u2(0)
u0=[1.5; 3];
% Lösningsintervall
T=8; tspan=linspace(0,T,200);
% Lös ekv.systemet med ode45 på interval tspan med begynnelsevärden u0.
[t,U]=ode45(f,tspan,u0)
% Tidsfunktion
subplot(1,3,1)
plot(t,U(:,1),t,U(:,2)), grid on, hold on, legend('u_1', 'u_2'), xlabel('t')
% Riktiningsfält
% First, set up graphics
u1=linspace(0,3,30);u2=linspace(0,3,30);
[U1,U2]=meshgrid(u1,u2);
% Then use quiver with the individual functions, and our grid as points, as
% parameters.
subplot(1,3,2:3)
quiver(U1,U2,f1(U1,U2),f2(U1,U2), 0.9), hold on
axis([0 3 0 3]), axis square, xlabel('u_1'), ylabel('u_2');

% Riktningsbanor
for	n1 = 0:0.5:3
    for n2 = 0:0.5:3
        u0(1)=n1;
        u0(2)=n2;
        [t,U]=ode45(f,tspan,u0)
        plot(U(:,1),U(:,2), 'r', 'LineWidth', 1.5)
    end
end
