```/* nag_ode_bvp_ps_lin_solve (d02uec) Example Program.
*
* Copyright 2019 Numerical Algorithms Group.
*
* Mark 27.0, 2019.
*/
#include <math.h>
#include <nag.h>

#ifdef __cplusplus
extern "C"
{
#endif
static double NAG_CALL exact(double x, Integer q);
static void NAG_CALL bndary(Integer m, double a, double b, double y[],
double bmat[], double bvec[]);
static void NAG_CALL pdedef(Integer m, double f[]);
#ifdef __cplusplus
}
#endif

int main(void)
{
/* Scalars */
Integer exit_status = 0;
double a = -nag_math_pi / 2.0, b = nag_math_pi / 2.0, resid;
Integer i, j, n, m = 3;
double teneps = 10.0 * nag_machine_precision;
/* Arrays */
double *bmat = 0, *bvec = 0, *f = 0, *uerr = 0, *y = 0, *c = 0, *f0 = 0,
*u = 0, *uc = 0, *x = 0;
/* NAG types */
Nag_Boolean reqerr = Nag_FALSE;
NagError fail;

INIT_FAIL(fail);

printf("nag_ode_bvp_ps_lin_solve (d02uec) Example Program Results\n\n");

/* Skip heading in data file */
scanf("%*[^\n] ");
scanf("%" NAG_IFMT "%*[^\n] ", &n);
if (!(u = NAG_ALLOC((n + 1) * (m + 1), double)) ||
!(f0 = NAG_ALLOC((n + 1), double)) ||
!(c = NAG_ALLOC((n + 1), double)) ||
!(uc = NAG_ALLOC((n + 1) * (m + 1), double)) ||
!(x = NAG_ALLOC((n + 1), double)) ||
!(bmat = NAG_ALLOC(m * (m + 1), double)) ||
!(bvec = NAG_ALLOC(m, double)) ||
!(f = NAG_ALLOC((m + 1), double)) ||
!(uerr = NAG_ALLOC((m + 1), double)) || !(y = NAG_ALLOC(m, double))
)
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}

/* Set up domain, boundary conditions and definition. */
bndary(m, a, b, y, bmat, bvec);
pdedef(m, f);

/* nag_ode_bvp_ps_lin_cgl_grid (d02ucc).
* Generate Chebyshev Gauss-Lobatto solution grid.
*/
nag_ode_bvp_ps_lin_cgl_grid(n, a, b, x, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_ode_bvp_ps_lin_cgl_grid (d02ucc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}

/* Set up problem right-hand sides for grid and transform. */
for (i = 0; i < n + 1; i++)
f0[i] = 2.0 * sin(x[i]) - 2.0 * cos(x[i]);

/* nag_ode_bvp_ps_lin_coeffs (d02uac).
* Coefficients of Chebyshev interpolating polynomial from function values f0
* on Chebyshev grid.
*/
nag_ode_bvp_ps_lin_coeffs(n, f0, c, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_ode_bvp_ps_lin_coeffs (d02uac).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}

/* nag_ode_bvp_ps_lin_solve (d02uec).
* Solve given boundary value problem on Chebyshev grid, in coefficient space
* using an integral formulation of the pseudospectral method.
*/
nag_ode_bvp_ps_lin_solve(n, a, b, m, c, bmat, y, bvec, f, uc, &resid,
&fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_ode_bvp_ps_lin_solve (d02uec).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}

/* nag_ode_bvp_ps_lin_cgl_vals (d02ubc).
* Obtain function values from coefficients of Chebyshev polynomial.
* Also obtain first- to third-derivative values.
*/
for (i = 0; i < m + 1; i++) {
nag_ode_bvp_ps_lin_cgl_vals(n, a, b, i, &uc[(n + 1) * i], &u[(n + 1) * i],
&fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_ode_bvp_ps_lin_cgl_vals (d02ubc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
}
/* Print solution. */
printf("Numerical Solution U and its first three derivatives\n\n");
printf("%7s%12s%12s%11s%11s\n", "x", "U", "Ux", "Uxx", "Uxxx");
for (i = 0; i < n + 1; i++)
printf("%10.4f %10.4f %10.4f %10.4f %10.4f\n", x[i], u[i], u[(n + 1) + i],
u[(n + 1) * 2 + i], u[(n + 1) * 3 + i]);

if (reqerr) {
for (i = 0; i < m + 1; i++)
uerr[i] = 0.0;

for (i = 0; i < n + 1; i++)
for (j = 0; j <= m; j++)
uerr[j] = MAX(uerr[j], fabs(u[(n + 1) * j + i] - exact(x[i], j)));

for (i = 0; i <= m; i++) {
printf("Error in the order %1" NAG_IFMT " derivative of U is < %8"
NAG_IFMT " * machine precision.\n", i,
10 * ((Integer) (uerr[i] / teneps) + 1));
}
}
END:
NAG_FREE(c);
NAG_FREE(f0);
NAG_FREE(u);
NAG_FREE(uc);
NAG_FREE(x);
NAG_FREE(bmat);
NAG_FREE(bvec);
NAG_FREE(f);
NAG_FREE(uerr);
NAG_FREE(y);
return exit_status;
}

static double NAG_CALL exact(double x, Integer q)
{
switch (q) {
case 0:
return cos(x);
break;
case 1:
return -sin(x);
break;
case 2:
return -cos(x);
break;
case 3:
return sin(x);
break;
}
return 0.0;
}

static void NAG_CALL bndary(Integer m, double a, double b, double y[],
double bmat[], double bvec[])
{
Integer i;
/* Boundary condition on left side of domain. */
for (i = 0; i < 2; i++)
y[i] = a;

y[2] = b;
/* Set up Dirichlet condition using exact solution at x = a. */
for (i = 0; i < m * (m + 1); i++)
bmat[i] = 0.0;
for (i = 0; i < 3; i++)
bmat[i] = 1.0;
for (i = 1; i < 3; i++)
bmat[m + i] = 2.0;
for (i = 1; i < 3; i++)
bmat[m * 2 + i] = 3.0;

bvec[0] = 0.0;
bvec[1] = 2.0;
bvec[2] = -2.0;
}

static void NAG_CALL pdedef(Integer m, double f[])
{
f[0] = 1.0;
f[1] = 2.0;
f[2] = 3.0;
f[3] = 4.0;
}
```