function [R]=Riccati(x,y0)

// Riccati: Solution of y'=x^2+y^2, x>=0
//                      y(0)=y0

C=y0*gamma(1/4)/gamma(3/4)/2;
x2=x.^2./2;
// R=x.*(C*besselj(-3/4,x2)+besselj(3/4,x2))./(besselj(-1/4,x2)-C*besselj(1/4,x2));

// Problem: Bessel functions of negative order not implemented in SCILAB
// therefore: calculation of besselj(-3/4,x2) and  besselj(-1/4,x2)
// via recurrence relation:
bess34=0.5*besselj(1/4,x2)./x2-besselj(5/4,x2);
bess14=1.5*besselj(3/4,x2)./x2-besselj(7/4,x2);
R=x.*(C*bess34+besselj(3/4,x2))./(bess14-C*besselj(1/4,x2));
