View difference between Paste ID: f741f87c9 and

-
+#include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
 #include <math.h>
 
 // overwrite u on success
 double* stepEuler(double* u, int const n, double const mu) {
         int i;
         double *un=malloc(sizeof(double)*n);
         if (!un) {
                 return 0;
         }
         if (!u || n < 3) {
                 free(un);
                 return 0;
         }
         for (i=1; i<(n-1); ++i) { // inner points
                 un[i]=u[i]+mu*(u[i-1]+u[i+1]-2*u[i]);
         }
         un[0]=u[0]+mu*2*(u[1]-u[0]); // zero-flux on boundaries
         un[n-1]=u[n-1]+mu*2*(u[n-2]-u[n-1]);
         memcpy(u,un,sizeof(double)*n);
         free(un);
         return u;
 }
 
 double* makeu0(int const n) {
         double *u0=0;
         u0=malloc(sizeof(double) * n);
         if (u0) {
                 int i;
                 for (i=0; i<n; ++i) {
                         double xi=M_PI*i/n;
                         u0[i]=sin(xi)*sin(xi);
                 }
         }
         return u0;
 }
 
 double sum_u(double const *u, int const n) {
         if ((u) && (n > 1)) {
                 double sum=0.0;
                 int i;
                 for (i=0; i<n; ++i) {
                         sum += u[i];
                 }
                 return sum;
         } else {
                 return 0.0;
         }
 }
 
 void print_u(double const *u, int const n) {
         if ((u) && (n > 1)) {
                 int i;
                 for (i=0; i<n; ++i) {
                         printf("%g,",u[i]);
                 }
                 printf("\n");
         }
 }
 
 int main(int argc, char *argv[]) {
         int n;
         double *u0, *ret;
         if (argc != 2) {
                 printf("usage: ceuler n>=3\n");
                 return 1;
         } else {
                 n=atoi(argv[1]);
                 if (n < 3) {
                         return 0;
                 }
         }
         u0=makeu0(n);
         ret=stepEuler(u0,n,0.5);
         if (ret) {
                 printf("%g\n",sum_u(u0,n));
                 free(u0);
         } else {
                 printf("failure\n");
                 free(u0);
                 return 2;
         }
 }