diff --git a/scripttest/Makefile b/scripttest/Makefile
new file mode 100644
index 0000000..9641619
--- /dev/null
+++ b/scripttest/Makefile
@@ -0,0 +1,4 @@
+integral: integral.c
+	mpicc integral.c -o integral -lm
+clean:
+	rm integral
diff --git a/scripttest/integral.c b/scripttest/integral.c
new file mode 100644
index 0000000..400f2af
--- /dev/null
+++ b/scripttest/integral.c
@@ -0,0 +1,118 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <stdbool.h>
+#include <mpi.h>
+
+const bool debug = false;
+const double a = 0.0;
+const double b = 1.0;
+
+/* Function declarations */
+void readinput(int argc, char* argv[], int iam, double* n_p);
+double trapezoidal_rule(double left_endpt, double right_endpt, int trap_count,
+   double base_len);
+double f(double x);
+
+int main(int argc, char** argv) {
+   int iam, np, local_n;
+   double n, h, local_a, local_b;
+   double local_int, total_int;
+   double start, finish, loc_elapsed, elapsed;
+
+   /* Star the MPI environment and discover who we are */
+   MPI_Init(NULL,NULL);
+   MPI_Comm_rank(MPI_COMM_WORLD, &iam);
+   MPI_Comm_size(MPI_COMM_WORLD, &np);
+
+   readinput(argc, argv, iam, &n); // Read the user input from the command line
+
+   /*Note: h and local_n are the same for all processes (load balance!)*/
+   h = (b-a)/n;     /* basis of each trapezoid */
+   local_n = n/np;  /* number of trapezoids per process */
+
+   /* Length of each process' interval of integration = local_n*h. */
+   local_a = a + iam*local_n*h;
+   local_b = local_a + local_n*h;
+   if (debug) printf("Process %d a = %f b = %f\n",iam,local_a,local_b);
+   /* Start a timer */
+   MPI_Barrier(MPI_COMM_WORLD);
+   start = MPI_Wtime();
+   /* We compute the local integral on each process using the local endpoints*/
+   local_int = trapezoidal_rule(local_a, local_b, local_n, h);
+   if (debug) printf("Process %d local integral = %f\n",iam,local_int);
+   finish = MPI_Wtime();
+   loc_elapsed = finish-start;
+   /* Do a max reduce to get the time */
+   MPI_Reduce(&loc_elapsed, &elapsed, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);
+
+   /* Do a sum reduce of the integrals computed on each process */
+   MPI_Reduce(&local_int, &total_int, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+
+
+   if (iam == 0) {
+      printf("Employing n = %.0f trapezoids\n", n);
+      printf("Trapezoid per process n = %d trapezoids\n", local_n);
+      printf("Number of processes: %d\n",np);
+      printf("of the integral of f(x) from %.0f to %.0f = %1.16f\n", a, b, total_int);
+      printf("Total elapsed time = %f milliseconds \n", elapsed * 1000);
+   }
+
+   /* Shut down MPI */
+   MPI_Finalize();
+
+   return 0;
+}
+
+void readinput(int argc, char* argv[], int iam, double* n_p){
+/*READINPUT gets the user input: the number of trapezoids
+ * Input args:   iam:  process rank in MPI_COMM_WORLD
+ *               np:  number of processes in MPI_COMM_WORLD
+ * Output args:  n_p:  pointer to number of trapezoids
+ */
+	if (iam == 0) {
+		if (argc!= 2){
+		    fprintf(stderr, "usage: srun %s <number of trapezoids> \n", argv[0]);
+            fflush(stderr);
+            *n_p = -1;
+		} else {
+			*n_p = atoi(argv[1]);
+		}
+	}
+	/* Broadcasts value of n to each process */
+	MPI_Bcast(n_p, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+	
+	// negative n ends the program
+    if (*n_p <= 0) {
+        MPI_Finalize();
+        exit(-1); // we return an error!
+    }
+}
+
+double trapezoidal_rule(double left_endpt, double right_endpt, int trap_count, double base_len) {
+/* TRAPEZOIDAL_RULE the serial function applying the trapezoidal rule
+ * Input:   left_endpt
+ *          right_endpt
+ *          trap_count
+ *          base_len
+ * Output:  Approximated value of the integral   
+ */
+   double estimate, x;
+   int i;
+
+   estimate = (f(left_endpt) + f(right_endpt))/2.0;
+   for (i = 1; i <= trap_count-1; i++) {
+      x = left_endpt + i*base_len;
+      estimate += f(x);
+   }
+   estimate = estimate*base_len;
+
+   return estimate;
+}
+
+
+double f(double x) {
+/*F function to be integrated */
+   return (4.0*sqrt(1.0 - x*x));
+}