diff --git a/.vscode/launch.json b/.vscode/launch.json
index e840c48..cee83f0 100644
--- a/.vscode/launch.json
+++ b/.vscode/launch.json
@@ -18,6 +18,14 @@
           "kind": "bin"
         }
       },
+      "initCommands": [
+        "command script import \"/home/alex/.rustup/toolchains/nightly-x86_64-unknown-linux-gnu/lib/rustlib/etc/lldb_rust_formatters.py\"",
+        "type summary add --no-value --python-function lldb_rust_formatters.print_val -x \".*\" --category Rust",
+        "type category enable Rust"
+      ],
+      "sourceLanguages": [
+        "rust"
+      ],
       "args": [],
       "cwd": "${workspaceFolder}"
     },
@@ -36,7 +44,15 @@
           "kind": "bin"
         }
       },
+      "initCommands": [
+        "command script import \"/home/alex/.rustup/toolchains/nightly-x86_64-unknown-linux-gnu/lib/rustlib/etc/lldb_rust_formatters.py\"",
+        "type summary add --no-value --python-function lldb_rust_formatters.print_val -x \".*\" --category Rust",
+        "type category enable Rust"
+      ],
       "args": [],
+      "sourceLanguages": [
+        "rust"
+      ],
       "cwd": "${workspaceFolder}"
     }
   ]
diff --git a/src/entity.rs b/src/entity.rs
index 6c93e02..08ce862 100644
--- a/src/entity.rs
+++ b/src/entity.rs
@@ -42,6 +42,7 @@ type PointVar = Var<Point2, Region2>;
 
 #[derive(Debug)]
 pub struct Point {
+    // pub id: i64,
     pub pos: PointVar,
 }
 
@@ -58,4 +59,8 @@ struct Line {
     p2: PointRef,
     len: ScalarVar,
     dir: ScalarVar,
-}
\ No newline at end of file
+}
+
+// struct System {
+//     points: Vec<Point>,
+// }
diff --git a/src/main.rs b/src/main.rs
index 8a0bdbf..28fd943 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -26,15 +26,15 @@ fn main() {
         );
     };
     print_points(&points);
-    let c1 = relation::Coincident {
-        p1: origin.clone(),
-        p2: p1.clone(),
-    };
+    // let c1 = relation::Coincident {
+    //     p1: origin.clone(),
+    //     p2: p1.clone(),
+    // };
     let c2 = relation::PointAngle::new_vertical(p1.clone(), p2.clone());
     let c3 = relation::PointAngle::new_horizontal(p2.clone(), p3.clone());
     let c4 = relation::AlignedDistance::new_vertical(p1.clone(), p2.clone(), 12.);
     let mut relations: Vec<Box<dyn Relation>> =
-        vec![Box::new(c1), Box::new(c2), Box::new(c3), Box::new(c4)];
+        vec![/*Box::new(c1),*/ Box::new(c2), Box::new(c3), Box::new(c4)];
     let mut constrained: Vec<Box<dyn Relation>> = Vec::new();
     let mut any_underconstrained = true;
     let mut any_constrained = true;
diff --git a/src/math.rs b/src/math.rs
index 0b7dd6e..9b96c6c 100644
--- a/src/math.rs
+++ b/src/math.rs
@@ -213,3 +213,391 @@ impl Region2 {
         }
     }
 }
+
+mod solve {
+
+    use std::collections::BTreeSet;
+    use std::fmt;
+    use std::iter::FromIterator;
+
+    use crate::math::Scalar;
+
+    // an unknown variable with an id
+    #[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord)]
+    struct Unknown(i64);
+
+    type UnknownSet = BTreeSet<Unknown>;
+
+    trait Unknowns {
+        fn unknowns(&self) -> UnknownSet;
+        fn has_unknowns(&self) -> bool;
+        fn has_unknown(&self, u: Unknown) -> bool;
+    }
+
+    impl Unknowns for Scalar {
+        fn unknowns(&self) -> UnknownSet {
+            UnknownSet::new()
+        }
+        fn has_unknowns(&self) -> bool {
+            false
+        }
+        fn has_unknown(&self, _: Unknown) -> bool {
+            false
+        }
+    }
+
+    impl Unknowns for Unknown {
+        fn unknowns(&self) -> UnknownSet {
+            FromIterator::from_iter(Some(*self))
+        }
+        fn has_unknowns(&self) -> bool {
+            true
+        }
+        fn has_unknown(&self, u: Unknown) -> bool {
+            *self == u
+        }
+    }
+
+    impl fmt::Display for Unknown {
+        fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
+            write!(f, "u{}", self.0)
+        }
+    }
+
+    #[derive(Clone, Debug, PartialEq)]
+    enum Expr {
+        Unkn(Unknown),
+        Const(Scalar),
+        Plus(Box<Expr>, Box<Expr>),
+        Neg(Box<Expr>),
+        Times(Box<Expr>, Box<Expr>),
+        Inv(Box<Expr>),
+    }
+
+    impl Unknowns for Expr {
+        fn unknowns(&self) -> UnknownSet {
+            use Expr::*;
+            match self {
+                Unkn(u) => u.unknowns(),
+                Const(_) => UnknownSet::default(),
+                Plus(l, r) | Times(l, r) => l.unknowns().union(&r.unknowns()).cloned().collect(),
+                Neg(e) | Inv(e) => e.unknowns(),
+            }
+        }
+        fn has_unknowns(&self) -> bool {
+            use Expr::*;
+            match self {
+                Unkn(u) => u.has_unknowns(),
+                Const(_) => false,
+                Plus(l, r) | Times(l, r) => l.has_unknowns() || r.has_unknowns(),
+                Neg(e) | Inv(e) => e.has_unknowns(),
+            }
+        }
+        fn has_unknown(&self, u: Unknown) -> bool {
+            use Expr::*;
+            match self {
+                Unkn(u1) => u1.has_unknown(u),
+                Const(_) => false,
+                Plus(l, r) | Times(l, r) => l.has_unknown(u) || r.has_unknown(u),
+                Neg(e) | Inv(e) => e.has_unknown(u),
+            }
+        }
+    }
+
+    impl Expr {
+        fn new_plus(e1: Expr, e2: Expr) -> Expr {
+            Expr::Plus(Box::new(e1), Box::new(e2))
+        }
+        fn new_times(e1: Expr, e2: Expr) -> Expr {
+            Expr::Times(Box::new(e1), Box::new(e2))
+        }
+        fn new_neg(e1: Expr) -> Expr {
+            Expr::Neg(Box::new(e1))
+        }
+        fn new_inv(e1: Expr) -> Expr {
+            Expr::Inv(Box::new(e1))
+        }
+        fn new_minus(e1: Expr, e2: Expr) -> Expr {
+            Expr::Plus(Box::new(e1), Box::new(Expr::new_neg(e2)))
+        }
+        fn new_div(e1: Expr, e2: Expr) -> Expr {
+            Expr::Times(Box::new(e1), Box::new(Expr::new_inv(e2)))
+        }
+
+        fn is_zero(&self) -> bool {
+            use Expr::*;
+            match self {
+                Const(c) => relative_eq!(*c, 0.),
+                _ => false,
+            }
+        }
+
+        fn is_one(&self) -> bool {
+            use Expr::*;
+            match self {
+                Const(c) => relative_eq!(*c, 1.),
+                _ => false,
+            }
+        }
+
+        fn simplify(self) -> Expr {
+            use Expr::*;
+            match self {
+                Plus(l, r) => match (l.simplify(), r.simplify()) {
+                    (Const(lc), Const(rc)) => Const(lc + rc),
+                    (Const(c), ref o) | (ref o, Const(c)) if relative_eq!(c, 0.) => o.clone(),
+                    (Times(l1, l2), Times(r1, r2)) => {
+                        if l2 == r2 {
+                            Expr::new_times(Expr::Plus(l1, r1), *l2).simplify()
+                        } else if l1 == r1 {
+                            Expr::new_times(Expr::Plus(l2, r2), *l1).simplify()
+                        } else if l1 == r2 {
+                            Expr::new_times(Expr::Plus(l2, r1), *l1).simplify()
+                        } else if l2 == r1 {
+                            Expr::new_times(Expr::Plus(l1, r2), *l2).simplify()
+                        } else {
+                            Expr::new_plus(Times(l1, l2), Times(r1, r2))
+                        }
+                    }
+                    (l, r) => Self::new_plus(l, r),
+                },
+                Times(l, r) => match (l.simplify(), r.simplify()) {
+                    (Const(lc), Const(rc)) => Const(lc * rc),
+                    (Const(c), ref o) | (ref o, Const(c)) if relative_eq!(c, 1.) => o.clone(),
+                    (Inv(ref den), ref num) | (ref num, Inv(ref den)) if *num == **den => Const(1.),
+                    (l, r) => Self::new_times(l, r),
+                },
+                Neg(v) => match v.simplify() {
+                    Const(c) => Const(-c),
+                    Neg(v) => *v,
+                    e => Self::new_times(Const(-1.), e),
+                },
+                Inv(v) => match v.simplify() {
+                    Const(c) => Const(1. / c),
+                    Inv(v) => *v,
+                    e => Self::new_inv(e),
+                },
+                e => e,
+            }
+        }
+
+        fn distrubute(self) -> Expr {
+            use Expr::*;
+            match self {
+                Plus(l, r) => Expr::new_plus(l.distrubute(), r.distrubute()),
+                Times(l, r) => match (*l, *r) {
+                    (Plus(l, r), o) | (o, Plus(l, r)) => Expr::new_plus(
+                        Expr::Times(Box::new(o.clone()), l),
+                        Expr::Times(Box::new(o), r),
+                    )
+                    .distrubute(),
+                    (l, r) => Expr::new_times(l, r),
+                },
+                Neg(v) => match *v {
+                    Plus(l, r) => Expr::new_plus(Neg(l).distrubute(), Neg(r).distrubute()),
+                    Times(l, r) => Expr::new_times(Neg(l).distrubute(), *r),
+                    Neg(v) => v.distrubute(),
+                    Inv(v) => Expr::new_inv(Neg(v).distrubute()),
+                    e => Expr::new_neg(e),
+                },
+                Inv(v) => match *v {
+                    Plus(l, r) => Expr::new_plus(Inv(l).distrubute(), Inv(r).distrubute()),
+                    Times(l, r) => Expr::new_times(Inv(l).distrubute(), Inv(r).distrubute()),
+                    Inv(v) => v.distrubute(),
+                    e => Expr::new_inv(e),
+                },
+                e => e,
+            }
+        }
+
+        fn reduce(self, for_u: Unknown) -> Expr {
+            use Expr::*;
+            match self {
+                Plus(l, r) => match (l.reduce(for_u), r.reduce(for_u)) {
+                    (Const(lc), Const(rc)) => Const(lc + rc),
+                    (l, r) => Self::new_plus(l, r),
+                },
+                Times(l, r) => match (l.reduce(for_u), r.reduce(for_u)) {
+                    (Const(lc), Const(rc)) => Const(lc * rc),
+                    (l, r) => Self::new_times(l, r),
+                },
+                Neg(v) => match v.reduce(for_u) {
+                    Unkn(u) if u == for_u => Expr::new_times(Const(-1.), Unkn(u)),
+                    e => Self::new_neg(e),
+                },
+                Inv(v) => match v.reduce(for_u) {
+                    e => Self::new_inv(e),
+                },
+                e => e,
+            }
+        }
+    }
+
+    impl fmt::Display for Expr {
+        fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
+            use Expr::*;
+            match self {
+                Unkn(u) => write!(f, "{}", u),
+                Const(c) => write!(f, "{}", c),
+                Plus(l, r) => write!(f, "({}) + ({})", l, r),
+                Times(l, r) => write!(f, "({}) * ({})", l, r),
+                Neg(e) => write!(f, "-({})", e),
+                Inv(e) => write!(f, "1 / ({})", e),
+            }
+        }
+    }
+
+    #[derive(Clone, Debug, PartialEq)]
+    struct Eqn(Expr, Expr);
+
+    impl Unknowns for Eqn {
+        fn unknowns(&self) -> UnknownSet {
+            self.0
+                .unknowns()
+                .union(&self.1.unknowns())
+                .cloned()
+                .collect()
+        }
+        fn has_unknowns(&self) -> bool {
+            self.0.has_unknowns() || self.1.has_unknowns()
+        }
+        fn has_unknown(&self, u: Unknown) -> bool {
+            self.0.has_unknown(u) || self.1.has_unknown(u)
+        }
+    }
+
+    fn ord_by_unkn(a: Expr, b: Expr, u: Unknown) -> Option<(Expr, Expr)> {
+        if a.has_unknown(u) {
+            Some((a, b))
+        } else if b.has_unknown(u) {
+            Some((b, a))
+        } else {
+            None
+        }
+    }
+
+    impl Eqn {
+        fn solve(&self, for_u: Unknown) -> Option<Expr> {
+            use Expr::*;
+            if !self.has_unknown(for_u) {
+                return None;
+            }
+            let (l, r) = (self.0.clone().simplify(), self.1.clone().simplify());
+            let (mut l, mut r) = ord_by_unkn(l, r, for_u)?;
+            loop {
+                let (new_l, new_r): (Expr, Expr) = match l {
+                    Unkn(u) => return if u == for_u { Some(r.simplify()) } else { None },
+                    Plus(a, b) => {
+                        let (a, b) = ord_by_unkn(*a, *b, for_u)?;
+                        (a, Expr::new_minus(r, b))
+                    }
+                    Times(a, b) => {
+                        let (a, b) = ord_by_unkn(*a, *b, for_u)?;
+                        (a, Expr::new_div(r, b))
+                    }
+                    Neg(v) => (*v, Expr::new_neg(r)),
+                    Inv(v) => (*v, Expr::new_inv(r)),
+                    Const(_) => return None,
+                };
+                l = new_l;
+                r = new_r;
+            }
+        }
+    }
+
+    #[derive(Clone, Debug, PartialEq)]
+    struct Eqns(Vec<Eqn>);
+
+    impl Unknowns for Eqns {
+        fn unknowns(&self) -> UnknownSet {
+            self.0.iter().flat_map(|eqn: &Eqn| eqn.unknowns()).collect()
+        }
+        fn has_unknowns(&self) -> bool {
+            self.0.iter().any(|eqn: &Eqn| eqn.has_unknowns())
+        }
+        fn has_unknown(&self, u: Unknown) -> bool {
+            self.0.iter().any(|eqn: &Eqn| eqn.has_unknown(u))
+        }
+    }
+
+    #[cfg(test)]
+    mod tests {
+        use super::*;
+
+        #[test]
+        fn test_unknowns() {
+            let u1 = Unknown(1);
+            let u2 = Unknown(2);
+            let u3 = Unknown(3);
+            assert!(u1.has_unknowns());
+            assert!(u2.has_unknowns());
+            assert!(u1.has_unknown(u1));
+            assert!(!u1.has_unknown(u2));
+            assert!(u1.unknowns().contains(&u1));
+            assert!(!u2.unknowns().contains(&u1));
+
+            let e1 = Expr::new_minus(Expr::Unkn(u1), Expr::Unkn(u2));
+            assert!(e1.has_unknowns());
+            assert!(e1.has_unknown(u1));
+            assert!(e1.has_unknown(u2));
+            assert!(!e1.has_unknown(u3));
+            assert!(e1.unknowns().len() == 2);
+        }
+
+        fn const_expr(e: Expr) -> Option<Scalar> {
+            match e {
+                Expr::Const(c) => Some(c),
+                _ => None,
+            }
+        }
+
+        #[test]
+        fn test_solve() {
+            use Expr::*;
+            let u1 = Unknown(1);
+            let e1 = Unkn(u1);
+            let e2 = Const(1.);
+
+            let eqn = Eqn(e1.clone(), e2.clone());
+            assert_eq!(eqn.solve(u1), Some(Const(1.)));
+            let eqn = Eqn(e2.clone(), e1.clone());
+            assert_eq!(eqn.solve(u1), Some(Const(1.)));
+
+            let e3 = Expr::new_plus(Const(1.), Const(1.));
+            let eqn = Eqn(e1.clone(), e3.clone());
+            assert_eq!(eqn.solve(u1), Some(Const(2.)));
+            let e3 = Expr::new_minus(Const(1.), Const(1.));
+            let eqn = Eqn(e1.clone(), e3.clone());
+            assert_eq!(eqn.solve(u1), Some(Const(0.)));
+
+            let e1 = Expr::new_div(Const(2.), Expr::new_minus(Const(1.), Const(4.)));
+            let e2 = Expr::new_minus(Const(1.), Unkn(u1));
+            let eqn = Eqn(e1, e2);
+            let e = eqn.solve(u1).unwrap();
+            assert!(const_expr(e.clone()).is_some());
+            assert!(relative_eq!(const_expr(e.clone()).unwrap(), 5. / 3.));
+
+            let e1 = Expr::new_times(Const(2.), Expr::new_minus(Const(1.), Const(4.)));
+            let e2 = Expr::new_minus(Expr::new_times(Unkn(u1), Const(2.)), Unkn(u1));
+            println!(
+                "e1==e2: {}=={} => {}=={}",
+                e1,
+                e2,
+                e1.clone().simplify(),
+                e2.clone().simplify()
+            );
+            println!(
+                "e1==e2: {}=={} => {}=={}",
+                e1,
+                e2,
+                e1.clone().distrubute(),
+                e2.clone().distrubute()
+            );
+            let eqn = Eqn(e1, e2);
+            let e = eqn.solve(u1).unwrap();
+            assert!(const_expr(e.clone()).is_some());
+            assert!(relative_eq!(const_expr(e.clone()).unwrap(), -6.));
+        }
+
+    }
+}