clear all
close all
clc

t0=0;
tk=0.001;
n=1000;
h=(tk-t0)/n;
t=t0:h:tk;

R1=1000;
C1=3.2e-8;
R2=1000;
C2=3.2e-8;
RL=1000;

for sym_no=1:7

    fu1=@(t,u1,u2) ((1/C1)*((1/R1)*(fe(t,sym_no)-u1)-((1/R1)*(u1-u2))));
    fu2=@(t,u1,u2) ((1/C2)*(((1/R2)*(u1-u2)-(u2/RL))));

    % zwykly euler

    eu1=[];
    eu2=[];
    eil=[];
    eu1(1)=0;
    eu2(1)=0;
    eil(1)=0;

    for i=1:n
        eu1(i+1)=eu1(i)+h*fu1(t(i),eu1(i),eu2(i));
        eu2(i+1)=eu2(i)+h*fu2(t(i),eu1(i),eu2(i));
        eil(i+1)=eu2(i+1)/RL;
    end

    wplot(t, eu1, eu2, eil, sym_no, 1);

    % ulepszony euler

    tu1=[];
    tu2=[];
    til=[];
    tu1(1)=0;
    tu2(1)=0;
    til(1)=0

    for i=2:n+1
        tu1(i)=tu1(i-1)+h*fu1((t(i-1)+(h/2)),(tu1(i-1)+(h/2)*fu1(t(i-1),tu1(i-1),tu2(i-1))),(tu2(i-1)+(h/2)*fu2(t(i-1),tu1(i-1),tu2(i-1))));
        tu2(i)=tu2(i-1)+h*fu2((t(i-1)+(h/2)),(tu1(i-1)+(h/2)*fu1(t(i-1),tu1(i-1),tu2(i-1))),(tu2(i-1)+(h/2)*fu2(t(i-1),tu1(i-1),tu2(i-1))));
        til(i)=tu2(i)/RL;
    end

    wplot(t, tu1, tu2, til, sym_no, 2);
end


function [e]=fe(t,sym_no)
    switch sym_no
        case 1
            e=2;
        case 2
            e=sin(2*pi*50*t);
        case 3
            e=sin(2*pi*600*t);
        case 4
            e=sin(2*pi*1750*t);
        case 5
            e=sin(2*pi*12000*t);
        case 6
            e=sin(2*pi*21000*t);
        case 7
            if mod(1000*t,5)<(1/2)
                e=1;
            elseif mod(1000*t,5)>=(1/2)
                e=0;
            end

    end
end

function [title]=wtitle(sym_no, euler_type)
    switch sym_no
        case 1
            title = 'Wykres dla E=2';
        case 2
            title = 'Wykres dla E=sin(2*pi*f*t), dla f=50 Hz';
        case 3
            title = 'Wykres dla E=sin(2*pi*f*t), dla f=600 Hz';
        case 4
            title = 'Wykres dla E=sin(2*pi*f*t), dla f=1750 Hz';
        case 5
            title = 'Wykres dla E=sin(2*pi*f*t), dla f=12 kHz';
        case 6
            title = 'Wykres dla E=sin(2*pi*f*t), dla f=21 kHz';
        case 7
            title = 'Wykres dla E=1 przy t<T/2 oraz E=0 przy t>=T/2 dla okresu 0.5ms';
    end

    if euler_type==1
        title = strcat(title, newline + "metoda Eulera");
    elseif euler_type==2
        title = strcat(title, newline + "metoda ulepszone Eulera");
    end
end

function wplot(t, u1, u2, il, sym_no, euler_type)
    figure
    plot(t,[u1;u2;il]);
    title(wtitle(sym_no, euler_type));
    legend('U1','U2','IL');
    xlabel('t (ms)');
    ylabel('V');
    hold off
end