Untitled

% main.m 
clc
clear all
close all

a = 0;
b = 6;
h = 0.5;
y0 = 1;

N = (b - a) / h;
x = zeros(N+1, 1);
for n = 1:(N+1)
    x(n) = a + (n-1)*h;
end

y_ex = zeros(N+1, 1);
y = @(x) (3/2 * x.^2 + 3/4 * x.^4 + 1) .^ (1/3);
for n = 1:(N+1)
    y_ex(n) = y(x(n));
end

f = @(x, y) x * (1 + x.^2) / y.^2;
y_ap = Xapxiptvpc1(a, b, f, h, y0);

plot(x, y_ex, 'r');
hold on;
plot(x, y_ap, 'b');
xlabel('x'); ylabel('y(x)');
legend('y_{ex}', 'y_{ap}');
title('Do thi nghiem chinh xac va nghiem xap xi cua bai toan Cauchy voi h=1/2');



% Xapxiptvpc1.m
function [y] = Xapxiptvpc1(a, b, f, h, y0)
    N = (b - a) / h;
    
    x = zeros(N+1, 1);
    for n = 1:(N+1)
        x(n) = a + (n-1)*h;
    end
    
    y = zeros(N+1, 1);
    y(1) = y0;
    for n = 1:N
        y(n+1) = y(n) + h * f(x(n), y(n));
    end
end