function [A,b,c] = rk4(c2,c3)
%  function [A,b,c] = rk4(c2,c3)
%  
%  Calculates coefficient arrays of four-stage 4th-order
%  Runge-Kutta formula depending on quadrature nodes c2, c3.

     c   = [0; c2; c3; 1];
     V   = flipud(vander(c)');
     r   = 1./[1:4]';
     b   = V \ r;
     a32 = c3*(c2-c3)/c2/(-1+2*c2)/2;
     a42 = (c2-1)*(-4*c3^2+5*c3+c2-2)/         ...
	   (-4*c2+6*c2*c3+3-4*c3)/(c2-c3)/c2/2;
     a43 = (-1+2*c2)*(c2-1)*(-1+c3)/           ...
           c3/(-4*c2+6*c2*c3+3-4*c3)/(c2-c3);
     A   = [  0           0    0  0;
             c2           0    0  0;
	     c3-a32     a32    0  0;
	     1-a42-a43  a42  a43  0];