<![CDATA[
#include <numeric.h>
#include <stdio.h>
#include <chplot.h>
#define NVAR 4
#define POINTS 50
int main() {
double t0=1, tf=10, y0[NVAR];
void func(double t, double y[], double dydt[], void *param);
double t[POINTS], y[NVAR][POINTS];
int i, points;
class CPlot plot;
y0[0]=j0(t0);
y0[1]=j1(t0);
y0[2]=jn(2,t0);
y0[3]=jn(3,t0);
points = oderk(func, t0, tf, y0, NULL, t, y);
printf("points = %d\n", points);
printf("\n%8s %18s %15s\n","t","integral","bessj(2,t)");
for (i=0; i<points; i++)
printf("%10.4f %14.6f %14.6f\n", t[i],y[2][i],jn(2,t[i]));
printf("\nThe leftover space in t and y are padded with results at tf\n");
for (i=points; i<POINTS; i++)
printf("%10.4f %14.6f %14.6f\n", t[i],y[2][i],jn(2,t[i]));
plotxy(t, y, "solving ode using oderk()", "t", "yi", &plot);
plot.legend("j0", 0);
plot.legend("j1", 1);
plot.legend("j2", 2);
plot.legend("j3", 3);
plot.plotting();
}
void func(double t, double y[], double dydt[], void *param) {
dydt[0] = -y[1];
dydt[1]=y[0]-(1.0/t)*y[1];
dydt[2]=y[1]-(2.0/t)*y[2];
dydt[3]=y[2]-(3.0/t)*y[3];
}
]]>
<![CDATA[
points = 43
t integral bessj(2,t)
1.0000 0.114903 0.114903
1.1995 0.159225 0.159225
1.4208 0.212473 0.212473
1.6107 0.259628 0.259628
1.8230 0.311687 0.311687
2.0595 0.365944 0.365944
2.3217 0.417776 0.417776
2.5117 0.447689 0.447689
2.7100 0.470487 0.470487
2.9266 0.484255 0.484255
3.1546 0.485094 0.485094
3.3763 0.471942 0.471942
3.5876 0.446661 0.446661
3.7842 0.412467 0.412467
4.0366 0.354910 0.354910
4.2242 0.303546 0.303546
4.4402 0.237267 0.237267
4.6658 0.162285 0.162285
4.8797 0.088340 0.088340
5.0780 0.019679 0.019679
5.3301 -0.064486 -0.064486
5.5121 -0.120926 -0.120926
5.7491 -0.186326 -0.186326
5.9414 -0.230973 -0.230973
6.1517 -0.269627 -0.269627
6.3458 -0.294896 -0.294896
6.5806 -0.311270 -0.311270
6.7618 -0.313089 -0.313089
6.9592 -0.304516 -0.304516
7.2117 -0.278417 -0.278417
7.3951 -0.249841 -0.249841
7.6042 -0.208797 -0.208797
7.8362 -0.154879 -0.154879
8.0631 -0.096271 -0.096271
8.2728 -0.039407 -0.039407
8.5323 0.030949 0.030949
8.7792 0.094212 0.094212
9.0451 0.154300 0.154300
9.2558 0.193600 0.193600
9.4769 0.225200 0.225200
9.6799 0.244472 0.244472
9.9231 0.254481 0.254481
10.0000 0.254630 0.254630
The leftover space in t and y are padded with results at tf
10.0000 0.254630 0.254630
10.0000 0.254630 0.254630
10.0000 0.254630 0.254630
10.0000 0.254630 0.254630
10.0000 0.254630 0.254630
10.0000 0.254630 0.254630
10.0000 0.254630 0.254630
]]>