include c:\cxpl\codes;          \intrinsic 'code' declarations

func real Func(FN, X);          \Return F(X) for function number FN
int  FN;  real X;
[case FN of
  1:  return X*X*X;
  2:  return 1.0/X;
  3:  return X
other return 0.0;
];

func Integrate(A, B, FN, N);    \Display area under curve for function FN
real A, B;  int FN, N;          \limits A, B, and number of slices N
real DX, X, Area;               \delta X
int  I;
[DX:= (B-A)/float(N);
X:= A;  Area:= 0.0;             \rectangular left
for I:= 1 to N do
    [Area:= Area + Func(FN,X)*DX;  X:= X+DX];
RlOut(0, Area);
X:= A;  Area:= 0.0;             \rectangular right
for I:= 1 to N do
    [X:= X+DX;  Area:= Area + Func(FN,X)*DX];
RlOut(0, Area);
X:= A+DX/2.0;  Area:= 0.0;      \rectangular mid point
for I:= 1 to N do
    [Area:= Area + Func(FN,X)*DX;  X:= X+DX];
RlOut(0, Area);
X:= A;  Area:= 0.0;             \trapezium
for I:= 1 to N do
    [Area:= Area + (Func(FN,X)+Func(FN,X+DX))/2.0*DX;  X:= X+DX];
RlOut(0, Area);
X:= A;  Area:= 0.0;             \Simpson's rule
for I:= 1 to N do
    [Area:= Area +
        DX/6.0*(Func(FN,X) + 4.0*Func(FN,(X+X+DX)/2.0) + Func(FN,X+DX));
    X:= X+DX];
RlOut(0, Area);
CrLf(0);
];

[Format(9,6);
Integrate(0.0, 1.0, 1, 100);
Integrate(1.0, 100.0, 2, 1000);
Integrate(0.0, 5000.0, 3, 5_000_000);
Integrate(0.0, 6000.0, 3, 6_000_000);
]