function hh=rk4 (hamso,x,is_plot,new_plot)
    %
    if (nargin <3) || isempty(is_plot)
        is_plot=0;
    end;
    if (nargin > 3) && ~isempty(new_plot) && (new_plot==1)
        figure();
    end;
    % de so sanh
    if is_plot==1
        %pp_euler(hamso,x,is_plot);
        rk2(hamso,x,is_plot);
    end;

    h=diff(x);
    [nouse g]=feval(hamso,x,0);
	y0=g(1)
    %y0=1.03709

    k1=zeros(size(x));
    k2=k1;
    k3=k1;
    k4=k1;
    
    print_tieude(27);
    for ii=2:length(x)
       [y0(ii) nouse k1(ii-1) k2(ii-1) k3(ii-1) k4(ii-1)]=pp_rk4_1b(hamso,x(ii-1),y0(ii-1),h(ii-1));
       str_print = sprintf('| %4.2g | %3.1f | %10.5f | %10.5f | %10.5f | %10.5f | %10.5f |',x(ii-1),h(ii-1),k1(ii-1),k2(ii-1),k3(ii-1),k4(ii-1),y0(ii));
       disp(str_print)
    end;
    if is_plot
		plot(x,g,'k-');
		hold on;
		plot(x,y0,'go','LineWidth',3);
		plot(x,y0,'g-.','LineWidth',3);
		grid on;
    end;
    %if (nargin > 3) && is_plot3d
    %    ham3chieu(hamso,x);
    %end;
function [f, g, k1, k2, k3, k4]=pp_rk4_1b(hamso,x0,y0,h)
	%
    [k1 g]=feval(hamso,x0,y0);
    
    k2 = feval(hamso,x0+h/2,y0+h/2*k1);
    k3 = feval(hamso,x0+h/2,y0+h/2*k2);
    k4 = feval(hamso,x0+h,y0+h*k3);
    f = y0 + h/6*(k1+2*k2+2*k3+k4);
    
function print_tieude(max_len)
    minus_str(1,1:(max_len*3+7+7))='-';
    disp(minus_str);
    str=printstr('x0',6);
    str=['|' str '|' printstr('h',5) '|'];
    str=[str printstr('k1=f(x0,y0)',13) '|'];
    str=[str printstr('k2=f(x0+h/2,y0+h/2*k1)',max_len) '|'];
    str=[str printstr('k3=f(x0+h/2,y0+h/2*k2)',max_len) '|'];
    str=[str printstr('k4=f(x0+h,y0)',15) '|'];
    printstr('y1=y0+h*(k1+2*k2+2*k3+k4)/6',max_len)
    str=[str printstr('y1=y0+h*(k1+2*k2+2*k3+k4)/6',max_len) '|'];
    disp(str);
    disp(minus_str);
  