diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
new file mode 100644
index 0000000..ae57107
--- /dev/null
+++ b/.github/workflows/test.yml
@@ -0,0 +1,13 @@
+name: "test pde"
+on:
+  push:
+  pull_request:
+
+jobs:
+  test:
+    name: cargo test
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions-rust-lang/setup-rust-toolchain@v1
+      - run: cargo test --all-features
\ No newline at end of file
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..ea8c4bf
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1 @@
+/target
diff --git a/Cargo.lock b/Cargo.lock
new file mode 100644
index 0000000..99f1c6c
--- /dev/null
+++ b/Cargo.lock
@@ -0,0 +1,182 @@
+# This file is automatically @generated by Cargo.
+# It is not intended for manual editing.
+version = 3
+
+[[package]]
+name = "autocfg"
+version = "1.3.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0c4b4d0bd25bd0b74681c0ad21497610ce1b7c91b1022cd21c80c6fbdd9476b0"
+
+[[package]]
+name = "black_scholes"
+version = "0.10.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "458d71a093446ab944668310ce42cad5fa704b411d95824d7e3b5169d65c7122"
+dependencies = [
+ "nrfind",
+ "serde",
+ "special",
+]
+
+[[package]]
+name = "libm"
+version = "0.2.8"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "4ec2a862134d2a7d32d7983ddcdd1c4923530833c9f2ea1a44fc5fa473989058"
+
+[[package]]
+name = "nrfind"
+version = "1.0.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7305df9ac15a14f878ba0b398faa954cad45d909827aac99757ed547a8a37fa2"
+dependencies = [
+ "num",
+]
+
+[[package]]
+name = "num"
+version = "0.2.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b8536030f9fea7127f841b45bb6243b27255787fb4eb83958aa1ef9d2fdc0c36"
+dependencies = [
+ "num-bigint",
+ "num-complex",
+ "num-integer",
+ "num-iter",
+ "num-rational",
+ "num-traits",
+]
+
+[[package]]
+name = "num-bigint"
+version = "0.2.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "090c7f9998ee0ff65aa5b723e4009f7b217707f1fb5ea551329cc4d6231fb304"
+dependencies = [
+ "autocfg",
+ "num-integer",
+ "num-traits",
+]
+
+[[package]]
+name = "num-complex"
+version = "0.2.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b6b19411a9719e753aff12e5187b74d60d3dc449ec3f4dc21e3989c3f554bc95"
+dependencies = [
+ "autocfg",
+ "num-traits",
+]
+
+[[package]]
+name = "num-integer"
+version = "0.1.46"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7969661fd2958a5cb096e56c8e1ad0444ac2bbcd0061bd28660485a44879858f"
+dependencies = [
+ "num-traits",
+]
+
+[[package]]
+name = "num-iter"
+version = "0.1.45"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1429034a0490724d0075ebb2bc9e875d6503c3cf69e235a8941aa757d83ef5bf"
+dependencies = [
+ "autocfg",
+ "num-integer",
+ "num-traits",
+]
+
+[[package]]
+name = "num-rational"
+version = "0.2.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5c000134b5dbf44adc5cb772486d335293351644b801551abe8f75c84cfa4aef"
+dependencies = [
+ "autocfg",
+ "num-bigint",
+ "num-integer",
+ "num-traits",
+]
+
+[[package]]
+name = "num-traits"
+version = "0.2.19"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841"
+dependencies = [
+ "autocfg",
+]
+
+[[package]]
+name = "option_price_pde"
+version = "0.1.0"
+dependencies = [
+ "black_scholes",
+]
+
+[[package]]
+name = "proc-macro2"
+version = "1.0.86"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5e719e8df665df0d1c8fbfd238015744736151d4445ec0836b8e628aae103b77"
+dependencies = [
+ "unicode-ident",
+]
+
+[[package]]
+name = "quote"
+version = "1.0.36"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0fa76aaf39101c457836aec0ce2316dbdc3ab723cdda1c6bd4e6ad4208acaca7"
+dependencies = [
+ "proc-macro2",
+]
+
+[[package]]
+name = "serde"
+version = "1.0.207"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5665e14a49a4ea1b91029ba7d3bca9f299e1f7cfa194388ccc20f14743e784f2"
+dependencies = [
+ "serde_derive",
+]
+
+[[package]]
+name = "serde_derive"
+version = "1.0.207"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "6aea2634c86b0e8ef2cfdc0c340baede54ec27b1e46febd7f80dffb2aa44a00e"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn",
+]
+
+[[package]]
+name = "special"
+version = "0.10.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b89cf0d71ae639fdd8097350bfac415a41aabf1d5ddd356295fdc95f09760382"
+dependencies = [
+ "libm",
+]
+
+[[package]]
+name = "syn"
+version = "2.0.74"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1fceb41e3d546d0bd83421d3409b1460cc7444cd389341a4c880fe7a042cb3d7"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "unicode-ident",
+]
+
+[[package]]
+name = "unicode-ident"
+version = "1.0.12"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3354b9ac3fae1ff6755cb6db53683adb661634f67557942dea4facebec0fee4b"
diff --git a/Cargo.toml b/Cargo.toml
new file mode 100644
index 0000000..715c2a7
--- /dev/null
+++ b/Cargo.toml
@@ -0,0 +1,9 @@
+[package]
+name = "option_price_pde"
+version = "0.1.0"
+edition = "2021"
+
+[dependencies]
+
+[dev-dependencies]
+black_scholes = "0.10.2"
diff --git a/src/lib.rs b/src/lib.rs
new file mode 100644
index 0000000..d9d75d3
--- /dev/null
+++ b/src/lib.rs
@@ -0,0 +1,245 @@
+//v is updated in place
+//a is length n-1, b length n, c length n-1, v length n
+//overwrite is appropriate if the matrix coefficients change (eg, if they depend on x), since they have to be recreated anyway
+fn thomas_algorithm_overwrite(a: &[f64], b: &[f64], c: &mut [f64], v: &mut [f64]) {
+    let mut solution = vec![0.0; v.len()];
+    c[0] = c[0] / b[0];
+    v[0] = v[0] / b[0];
+
+    for i in 1..c.len() {
+        let den = b[i] - a[i - 1] * c[i - 1];
+        c[i] = c[i] / den;
+        v[i] = (v[i] - a[i - 1] * v[i - 1]) / den;
+    }
+    let last_index = v.len() - 1;
+    v[last_index] = (v[last_index] - a[last_index - 1] * v[last_index - 1])
+        / (b[last_index] - a[last_index - 1] * c[last_index - 1]);
+
+    solution[last_index] = v[last_index];
+
+    for i in (0..last_index).rev() {
+        v[i] = v[i] - c[i] * v[i + 1];
+    }
+}
+
+pub struct Surface {
+    pub surface: Vec<Vec<f64>>,
+    pub size_t: usize,
+    pub size_x: usize,
+    pub dx: f64,
+    pub dt: f64,
+    pub min_t: f64,
+    pub min_x: f64,
+}
+impl Surface {
+    pub fn init(size_t: usize, size_x: usize, dt: f64, dx: f64, min_t: f64, min_x: f64) -> Self {
+        Surface {
+            surface: vec![],
+            size_t,
+            size_x,
+            dx,
+            dt,
+            min_t,
+            min_x,
+        }
+    }
+    //note that this forces a move
+    pub fn set_next_x(&mut self, x_at_t: Vec<f64>) {
+        self.surface.push(x_at_t);
+    }
+    pub fn get_last_x(&self) -> Option<&[f64]> {
+        self.get_x_at_t(self.surface.len() - 1)
+    }
+    pub fn get_x_at_t(&self, t_index: usize) -> Option<&[f64]> {
+        let curr_length = self.surface.len();
+        if t_index < curr_length {
+            Some(&self.surface[t_index])
+        } else {
+            None
+        }
+    }
+    pub fn get_result_at_x_and_t(&self, t_index: usize, x_index: usize) -> Option<(f64, f64, f64)> {
+        let x = self.get_x_at_t(t_index)?;
+        let curr_length = x.len();
+        if x_index < curr_length {
+            Some((
+                self.min_t + self.dt * (t_index as f64 + 1.0),
+                self.min_x + self.dx * (x_index as f64 + 1.0),
+                x[x_index],
+            ))
+        } else {
+            None
+        }
+    }
+}
+pub struct Boundary {
+    min: f64,
+    max: f64,
+}
+pub fn solve_second_order_pde(
+    fn_fn: impl Fn(f64, f64) -> f64, //for black scholes, this is a constant (-r)
+    fn_first_deriv: impl Fn(f64, f64) -> f64, //for black scholes, this is r*S (or, in logs, a constant r-sigma*sigma/2)
+    fn_second_deriv: impl Fn(f64, f64) -> f64, //for black scholes, this is sigma*sigma*S*S/2
+    boundary_up: impl Fn(f64, f64) -> f64,    //for black scholes, this is s-ke^{-rt}
+    boundary_down: impl Fn(f64, f64) -> f64,  //for black scholes, this is 0
+    init_fn: impl Fn(f64) -> f64,             //for black scholes, this is (s-k)^+
+    size_t: usize,
+    t_lim: &Boundary,
+    size_x: usize,
+    x_lim: &Boundary,
+) -> Option<Surface> {
+    let dx = (x_lim.max - x_lim.min) / (size_x as f64 + 1.0); //+ two so that boundary is "outside" of mesh
+    let dt = (t_lim.max - t_lim.min) / (size_t as f64); //+ one so that boundary is "outside" of mesh
+    let mut surface = Surface::init(size_t, size_x, dt, dx, t_lim.min, x_lim.min);
+    let mut construct_starting_surface = vec![0.0; size_x]; //does not include boundaries
+    for i in 0..size_x {
+        construct_starting_surface[i] = init_fn(x_lim.min + dx * (i as f64 + 1.0));
+        //start at 1 since 0 is boundary
+    }
+
+    let fxminus =
+        |t: f64, x: f64| (dt / dx) * (0.5 * fn_first_deriv(t, x) - fn_second_deriv(t, x) / dx);
+    let fxplus =
+        |t: f64, x: f64| -(dt / dx) * (0.5 * fn_first_deriv(t, x) + fn_second_deriv(t, x) / dx);
+    let fx =
+        |t: f64, x: f64| 1.0 + 2.0 * fn_second_deriv(t, x) * (dt / (dx * dx)) - fn_fn(t, x) * dt;
+
+    let put_boundarys = |v: &mut [f64], t: f64| {
+        v[0] = v[0] - fxminus(t, x_lim.min + dx) * boundary_down(t, x_lim.min); //a at x+1 for the "first" row in mesh
+        v[size_x - 1] = v[size_x - 1] - fxplus(t, x_lim.max - dx) * boundary_up(t, x_lim.max);
+        //c at x-1 for the "last" row in mesh
+    };
+    surface.set_next_x(construct_starting_surface);
+    //t=0 is the start of the surface, continuing with t=mint+dt
+    for j in 1..size_t {
+        let t = t_lim.min + dt * (j as f64);
+        let a: Vec<f64> =
+            (1..size_x) //start at 1 since 0 is boundary
+                .map(|index| fxminus(t, x_lim.min + dx * (index as f64 + 1.0)))
+                .collect(); //add 1 since the first "a" is actually the second row
+
+        let b: Vec<f64> = (1..size_x + 1) //start at 1 since 0 is boundary, end at size_x since size_x+1 is boundary
+            .map(|index| fx(t, x_lim.min + dx * (index as f64)))
+            .collect();
+        let mut c: Vec<f64> = (1..size_x) //start at 1 since 0 is boundary
+            .map(|index| fxplus(t, x_lim.min + dx * (index as f64)))
+            .collect();
+
+        let mut v = surface.get_last_x()?.to_vec().clone();
+        put_boundarys(&mut v, t);
+        thomas_algorithm_overwrite(&a, &b, &mut c, &mut v);
+        surface.set_next_x(v);
+    }
+    Some(surface)
+}
+
+#[cfg(test)]
+mod tests {
+
+    use super::*;
+    extern crate black_scholes;
+
+    #[test]
+    fn thomas_algorithm() {
+        let a = vec![2.0, -4.0];
+        let b = vec![1.0, 2.0, 1.0];
+        let mut c = vec![3.0, -1.0];
+        let mut d = vec![6.0, 4.0, 0.0];
+        /*let result = */
+        thomas_algorithm_overwrite(&a, &b, &mut c, &mut d);
+        for res in d.iter() {
+            println!("res {}", res);
+        }
+        assert_eq!(d[0], 3.0);
+        assert_eq!(d[1], 1.0);
+        assert_eq!(d[2], 4.0);
+    }
+
+    #[test]
+    fn test_pde_on_black_scholes() {
+        let stock = 5.0;
+        let strike = 4.5;
+        let rate = 0.01;
+        let sigma = 0.3;
+        let maturity = 2.0;
+        let t_lim = Boundary {
+            min: 0.0,
+            max: maturity,
+        };
+        let x_lim = Boundary {
+            min: -10.0,
+            max: 10.0,
+        };
+        let size_t = 700;
+        let size_x = 255;
+        let result = solve_second_order_pde(
+            |_t, _x| -rate,
+            |_t, _x| rate - sigma * sigma * 0.5,
+            |_t, _x: f64| sigma * sigma * 0.5,
+            |t, x| (stock * x.exp() - strike) * (-rate * t).exp(),
+            |_t, _x| 0.0,
+            |x| {
+                let v = stock * x.exp() - strike;
+                if v > 0.0 {
+                    v
+                } else {
+                    0.0
+                }
+            },
+            size_t,
+            &t_lim,
+            size_x,
+            &x_lim,
+        )
+        .unwrap();
+        assert_eq!(result.surface.len(), size_t); //starts at "0"
+        let (_t, x, approx_bs) = result
+            .get_result_at_x_and_t(size_t - 1, ((size_x - 1) / 2) as usize)
+            .unwrap();
+        assert!(
+            (approx_bs - black_scholes::call(stock * x.exp(), strike, rate, sigma, maturity)).abs()
+                < 0.01
+        );
+    }
+    #[test]
+    fn test_pde_on_full_range_bs() {
+        let strike = 4.5;
+        let rate = 0.01;
+        let sigma = 0.3;
+        let maturity = 2.0;
+        let t_lim: Boundary = Boundary {
+            min: 0.0,
+            max: maturity,
+        };
+        let x_lim = Boundary { min: 0.0, max: 1.0 };
+        let size_t = 700;
+        let size_x = 255;
+        let result = solve_second_order_pde(
+            |_t, x| -rate * (1.0 - x),
+            |_t, x| rate * (1.0 - x) * x,
+            |_t, x: f64| (1.0 - x) * (1.0 - x) * x * x * sigma * sigma * 0.5,
+            |_t, x| (2.0 * x - 1.0), // * (-rate * t).exp(),
+            |_t, _x| 0.0,
+            |x| {
+                let v = 2.0 * x - 1.0;
+                if v > 0.0 {
+                    v
+                } else {
+                    0.0
+                }
+            },
+            700,
+            &t_lim,
+            255,
+            &x_lim,
+        )
+        .unwrap();
+        let (_t, x, v) = result
+            .get_result_at_x_and_t(size_t - 1, ((size_x - 1) / 2) as usize)
+            .unwrap();
+        let s = (x * strike) / (1.0 - x);
+        let approx_bs = (s + strike) * v;
+        assert!((approx_bs - black_scholes::call(s, strike, rate, sigma, maturity)).abs() < 0.001);
+        //note that this is more accurate than the approximation of the boundary in the previous test
+    }
+}