% Contoh Program untuk mencari solusi persamaan differensial
% (yang ada di berkas dfungsi.m) dengan Metode Numerik GABUNGAN
% Error tidak dihitung dengan solusi analitik,
% Error di-estimasi dengan selisih solusi Order Pertama dan Kedua
%
To = input('Masukkan waktu awal To, biasanya 0 saja = ')
Ta = input('Masukkan waktu akhir Ta =')
TOL = input('Masukkan nilai toleransi (dalam %) = ')
N = input('Masukkan maksimum iterasi =')
delta_t = input('Masukkan nilai interval awal delta_t = ')
for I = 1:3 % diberi kesempatan 3 kali mencoba
    if delta_t > (Ta-To)/2
        'delta_t terlalu besar'
        delta_t = input('Masukkan lagi nilai interval awal delta_t = ')
    elseif delta_t < (Ta-To)/N
        'delta_t terlalu kecil'
        delta_t = input('Masukkan lagi nilai interval awal delta_t = ')
    end
end
Xo = input('masukkan kondisi awal x(0) = ')
K = round((Ta-To)/delta_t);
x(1) = Xo;
x1(1) = Xo;
t(1) = To;
% Tidak perlu dihitung: xANA(1) = x_analitik(t(1)); % Solusi analitik
galat(1) = 0; % Error dalam %
for i = 1:N
    t(i+1)= t(i) + delta_t;
    % Order Pertama
    x1(i+1) = x(i) + (delta_t*dfungsi(t(i),x(i)));
    % Order Kedua:
    % Dapat dihitung x'(i):
    dx(i) = dfungsi(t(i),x(i));
    % Dengan Metode Euler Order Pertama, dihitung estimasi Ex(i+1):
    Ex(i+1) = x(i) + (delta_t*dx(i));
    % Lalu dihitung estimasi x'(i+1):
    Edx(i+1) = dfungsi(t(i+1), Ex(i+1));
    x(i+1) = x(i)+0.5*delta_t*(dx(i) + Edx(i+1));
    % Estimasi Error:
    galat(i+1) = 100*abs((x1(i+1)-x(i+1))/x(i+1)); % Estimasi Error dalam %
    if galat(i+1) > TOL % Jika error masih lebih besar dari toleransi
        t(i+1) = t(i);
        delta_t = delta_t/10; % delta_t dibagi 10
    elseif delta_t <= (Ta-To)/N
        'Toleransi TOL terlalu ketat'
        DELTA_T = delta_t
        Iterasi_ke = i
        time = t(i+1)
        ERROR = galat(i+1)
        break
    elseif t(i+1) > Ta
        'Toleransi TOL terlalu longgar'
        DELTA_T = delta_t
        Iterasi_ke = i
        time = t(i+1)
        ERROR = galat(i+1)
        break
    end
end
subplot(211), plot(t,x,'go',t,x1,'r*'), grid on, ylabel('x_1 dan x_2')
subplot(212), plot(t,galat,'ro'), grid on, ylabel('ERROR'), xlabel('t')