function R=N85chapt9ex4Rungkuta4(a,b,N,ya)
%學號 20042374  姓名 鐘健升   作業號  85
%p264 第4題四階RK方法
%f是微分方程右端表達式
%g是微分方程解的表達式
%a,b是區間[a,b]中的a,b
%N是等分區間的個數
syms x y
f=1-2*x*y/(1+x^2);
g=(x+1/3*x^3)/(1+x^2);
h=(b-a)/N;
y=zeros(1,N+1);
x=zeros(1,N+1);
y(1)=ya;
x=a:h:b;
for i=1:N
    k1=1-2*x(i)*y(i)/(1+x(i)^2)
    k2=1-2*(x(i)+h/2)*(y(i)+(h/2)*k1)/(1+(x(i)+h/2)^2)
    k3=1-2*(x(i)+h/2)*(y(i)+(h/2)*k2)/(1+(x(i)+h/2)^2)
    k4=1-2*(x(i)+h)*(y(i)+h*k3)/(1+(x(i)+h)^2)
    y(i+1)=y(i)+(h/6)*(k1+2*k2+2*k3+k4);
end
E=[x',y']