Code:

// PROGRAM tdse
// motion of a wavepacket incident on a potential
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "area.h"
const double PI = 3.14159265358979323846264338327950;
static double Re[1101], Im[1101], Imold[1101];
void parameters(double *x0, double *k0, double *width, double *V0, double *a,
double *xmin, double *xmax, double *n, double *dx, double *dx2,
double *dt);
void initial_packet(double Re[], double Im[], double x0, double k0,
double width, double xmin, double n, double dx, double dt);
void set_up_windows(double xmin, double xmax, double V0,
AREA **area1, AREA **area2, AREA **area3, AREA **area4);
void draw_potential(double V0, double a, double xmin, double xmax,
AREA *area4);
int main(void)
{
void evolve(double Re[], double Im[], double Imold[], double *t,
double V0, double a, double dx, double dx2, double dt,
double xmin, double n);
void plots(double Re[], double Im[], double Imold[], double *t,
double xmin, double n, double dx, AREA *area1, AREA *area2,
AREA *area3, AREA *area4);
double x0, k0, width, V0, a, xmin, xmax, n, dx, dx2, dt;
parameters(&x0, &k0, &width, &V0, &a, &xmin, &xmax, &n, &dx, &dx2, &dt);
initial_packet(Re, Im, x0, k0, width, xmin, n, dx, dt);
AREA *area1, *area2, *area3, *area4;
set_up_windows(xmin, xmax, V0, &area1, &area2, &area3, &area4);
draw_potential(V0, a, xmin, xmax, area4);
double t = 0.0;
do {
evolve(Re, Im, Imold, &t, V0, a, dx, dx2, dt, xmin, n);
plots(Re, Im, Imold, &t, xmin, n, dx, area1, area2, area3, area4);
} while (getchar() != EOF);
closearea(area1);
closearea(area2);
closearea(area3);
closearea(area4);
return EXIT_SUCCESS;
}
void parameters(double *x0, double *k0, double *width, double *V0, double *a,
double *xmin, double *xmax, double *n, double *dx, double *dx2,
double *dt)
{
*x0 = -15;
*width = 1;
*k0 = 2;
*xmax = 20;
*xmin = -*xmax;
*V0 = 2;
*a = 1;
*dx = 0.4;
*dx2 = (*dx)*(*dx);
*n = (*xmax - *xmin)/(*dx);
*dt = 0.1;
}
void initial_packet(double Re[], double Im[], double x0, double k0,
double width, double xmin, double n, double dx, double dt)
// initial Gaussian wavepacket
{
double delta2 = width*width;
double A = pow(2*PI*delta2, -0.25);
double b = 0.5*k0*dt;
for (int i=1; i <= n; i++) {
double x = xmin + (i-1)*dx;
double arg = 0.25*pow(x - x0, 2)/delta2;
double e = exp(-arg);
Re[i] = A*cos(k0*(x - x0))*e;
Im[i] = A*sin(k0*(x - x0 - 0.5*b))*e;
}
}
void set_up_windows(double xmin, double xmax, double V0,
AREA **area1, AREA **area2, AREA **area3, AREA **area4)
{
*area1 = openarea(0, 1, 0.75, 1.0);
setwindow(*area1, xmin, xmax, -0.1, 0.5);
plotline(*area1, xmin, 0, xmax, 0);
*area2 = openarea(0, 1, 0.5, 0.75);
setwindow(*area2, xmin, xmax, -1, 1);
plotline(*area2, xmin, 0, xmax, 0);
*area3 = openarea(0, 1, 0.25, 0.5);
setwindow(*area3, xmin, xmax, -1, 1);
plotline(*area3, xmin, 0, xmax, 0);
*area4 = openarea(0, 1, 0, 0.22);
setwindow(*area4, xmin, xmax, -V0, V0);
}
void draw_potential(double V0, double a, double xmin, double xmax,
AREA *area4)
{
double V(double x, double V0, double a);
AREA *area = area4;
setcolor(area, "red");
plotline(area, xmin, 0, a, 0);
plotline(area, a, 0, a, V0);
plotline(area, a, V0, xmax, V0);
setcolor(area, "black/white");
printf("total probability = ");
}
void plots(double Re[], double Im[], double Imold[], double *t, double xmin,
double n, double dx, AREA *area1, AREA *area2, AREA *area3, AREA *area4)
{
AREA *area = area2;
cleararea(area);
double xmax = 0.0;
plotline(area, xmin, 0, xmax, 0);
for (int i=0; i <= n; i++) {
double x = xmin + (i - 1)*dx;
plotpoint(area, x, Im[i]);
}
area = area1;
cleararea(area);
plotline(area, xmin, 0, xmax, 0);
double Psum = 0;
for (int i=1 ; i <= n; i++) {
double x = xmin + (i - 1)*dx;
double P = Re[i]*Re[i] + Im[i]*Imold[i];
Psum += P*dx;
plotpoint(area, x, P);
}
area = area4;
setcursor(area, 1, 20);
printf(" ");
setcursor(area, 1, 20);
printf("%g", Psum);
}
void evolve(double Re[], double Im[], double Imold[], double *t,
double V0, double a, double dx, double dx2, double dt,
double xmin, double n)
{
double V(double x, double V0, double a);
for (int i=1; i <= n; i++) {
double x = xmin + (i-1)*dx;
double HIm = V(x,V0,a)*Im[i] - 0.5*(Im[i+1] - 2*Im[i] + Im[i-1])
/ dx2;
// real part defined at multiples of dt
Re[i] += HIm*dt;
}
for (int i=1; i <= n; i++) {
double x = xmin + (i-1)*dx;
// dt/2 earlier than real part
Imold[i] = Im[i];
double HRe = V(x,V0,a)*Re[i] - 0.5*(Re[i+1] - 2*Re[i] + Re[i-1])
/ dx2;
// dt/2 later than real part
Im[i] -= HRe*dt;
}
*t += dt; // time of real part
}
double V(double x, double V0, double a)
{
// step potential
if (x > a) {
return V0;
}else{
return 0;
}
}

If you want the graphics in C, you'll have to implement a compatible set of calls that manipulate the AREA objects. I don't really know how they're supposed to work, so these are just stub functions to get the thing to compile.